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Abstract 


In  the  large-eddy  simulation  (LES)  of  high  Reynolds  number  wall-bounded  flows,  wall  modeling 
is  needed  to  alleviate  the  severe  near- wall  resolution  requirement.  Simple  algebraic  models  such  as 
the  instantaneous  log-law  are  inadequate  for  predicting  complex  flows  with  strong  pressure  gradients 
and  separation.  Under  the  sponsorship  of  AFOSR,  we  have  explored  two  classes  of  wall  models: 
those  based  on  the  turbulent  boundary-layer  (TBL)  equations  and  those  based  on  control  theory. 
In  this  report,  a  recent  application  of  the  TBL  equation  model  to  LES  of  the  flow  over  a  cylinder  at 
super-critical  Reynolds  number  will  first  be  discussed.  The  emphasis  of  the  report  is,  however,  on 
control  based  wall  modeling,  in  which  sub-optimal  control  strategy  is  used  to  find  the  wall  stresses 
that  will  force  the  outer  LES  toward  a  target  profile.  Results  from  channel-flow  simulation  indicate 
that  in  order  to  obtain  the  correct  mean  velocity  profile  (the  log  law),  the  wall  stresses  must  not 
only  model  the  physics  but  also  compensate  for  numerical  and  SGS  modeling  errors.  The  data 
generated  by  this  sub-optimal  control  strategy  are  then  used  to  derive  a  linear  stochastic  estimate 
model.  The  mathematical  formulation  and  issues  of  key  importance  in  control-based  wall  modeling 
will  be  discussed  in  detail.  Current  efforts  towards  a  predictive  and  inexpensive  wall  model  in  the 
control  framework  will  also  be  presented. 
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Chapter  1 

Introduction 


Large-eddy  simulations  (LES)  of  high  Reynolds  number  flows  are  difficult  to  perform  due  to  the 
need  to  include  a  large  number  of  grid  points  in  the  near  wall  region.  While  LES  models  the  small 
scales  of  the  flow  and  resolves  the  large,  dynamically  important  scales,  near  the  wall,  eddies  scale 
with  the  distance  from  the  wall  and  move  increasingly  closer  to  the  wall  as  the  Reynolds  number 
increases.  These  eddies  are  dynamically  important  despite  their  size.  Unfortunately,  the  eddy 
viscosity  sub-grid  scale  (SGS)  models  only  make  a  small  contribution  to  the  total  Reynolds  stress. 
This  makes  these  models  invalid  near  the  wall  [21],  unless  the  LES  grid  is  sufficiently  refined  to 
resolve  the  near- wall  vortical  structures.  The  required  number  of  grid  points  for  such  a  wall-resolved 
LES  scales  as  Re \  in  an  attached  boundary  layer  [4],  which  is  only  a  slight  improvement  on  the 
scaling  for  a  fall  direct  numerical  simulation  (DNS)  of  iie9/4. 

The  technique  of  wall  modeling  was  developed  to  reduce  the  Reynolds  number  scaling  of  LES 
resolution,  so  that  LES  could  be  applied  in  practical  situations.  For  recent  reviews,  see  [13]  and  [35]. 
The  approach  has  a  long  history  dating  back  to  the  atmospheric  science  and  oceanographic  appli¬ 
cations.  Limited  by  the  computational  power  of  the  time,  Deardorff  [16]  was  the  first  to  implement 
a  model  for  the  wall  layer  in  an  LES  of  a  channel  flow  at  infinite  Reynolds  number.  He  imposed 
constraints  on  wall-parallel  velocities  in  terms  of  their  wall-normal  second  derivatives  to  ensure 
the  LES  satisfied  the  log-law  in  mean.  The  wall  transpiration  velocity  was  set  to  zero.  The  first 
“modern”  wall  model  was  developed  by  Schumann  [37].  It  is  a  modern  wall  model  in  the  sense 
that  the  wall  stresses  were  determined  directly  by  an  algebraic  model.  The  wall  stresses  were  found 
by  assuming  that  they  were  in  phase  with  the  velocity  at  the  first  off-wall  grid  point  and  that 
the  deviation  from  their  mean  was  proportional  to  the  deviation  of  the  velocity  from  its  mean. 
Since  the  flow  was  in  a  channel,  both  the  mean  wall  stresses  and  mean  velocities  were  known.  The 
transpiration  velocity,  as  in  the  case  of  Deardorff  [16],  was  set  to  zero. 

Many  improvements  to  this  basic  model  have  been  made  in  the  intervening  years.  For  example, 
Piomelli  et  al.  [36]  added  an  offset  to  the  velocity  to  account  for  the  inclination  of  the  eddies  to 
obtain: 


r12(x>  *) 

u(x  +  As,yi,z) 

(Tw)  > 

(u) 

V(x,0,  z) 

=  o, 

Tyiixi z) 

w(x  +  A  s,yi,z) 

(rw)  > 

(u) 

where  As  is  an  empirical  displacement  related  to  the  mean  orientation  of  vortical  structures  near 
the  wall,  yi  is  the  wall-normal  coordinate  of  the  first  off-wall  grid  points,  and  (*)  denotes  averaging 
in  the  wall-parallel  ( x-z )  plane.  For  atmospheric  boundary  layers,  Mason  &  Callen  [27]  enforces  an 
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instantaneous  logarithmic  law  in  each  grid  cell  adjacent  to  the  wall.  Other  efforts  include  those  by 
Grotzbach  [18],  Werner  &  Wengle  [44],  and  Hoffman  &  Benocci  [20].  While  these  types  of  algebraic 
models  produced  reasonable  results  in  the  LES  of  plane  channels  and  annuli,  they  are  in  general 
inadequate  for  flows  in  more  complex  geometries.  In  these  flows,  adverse  pressure  gradients  and 
separated  regions  are  common,  and  methods  based  on  an  equilibrium  balance  of  stresses  and  a 
logarithmic  velocity  profile  give  inaccurate  results. 

To  address  this  robustness  issue  in  wall  modeling,  several  investigators  used  more  elaborate 
near- wall  flow  models  to  compute  the  wall  stresses  (see  e.g.  [8]  and  [13]).  This  type  of  approach 
divides  the  computational  domain  into  two  regions:  one  near  the  wall  and  one  away  from  the  wall. 
A  simplified  set  of  equations  based  on  turbulent  boundary-layer  (TBL)  approximations  are  solved 
on  a  near  wall  grid  separate  from  the  outer  LES  grid,  subject  to  boundary  conditions  determined 
from  the  outer  LES  velocity  together  with  the  no-slip  wall.  The  equations  solved  in  the  wall-layer 
are  of  the  following  general  form, 


d_ 

dy 


(  »  \ 


1  dp  dui  d 
pdxi  dt  dxjU%U^ 


i  =  M, 


(1.2) 


although  simplified  versions,  with  the  right  hand  side  set  to  zero  or  the  pressure  gradient  alone, 
have  also  been  considered  [13,  42].  In  Eq.  (1.2)  vt  is  given  by  a  RANS  eddy  viscosity  model,  and 
the  pressure  is  imposed  from  the  LES  solution  and  assumed  constant  across  the  wall  layer.  The 
computed  wall  stress  is  then  provided  to  the  LES  as  a  boundary  condition.  While  this  method 
does  require  the  solution  of  an  extra  set  of  equations,  the  simplifications  made  in  these  equations 
makes  its  cost  much  less  than  the  evaluation  of  the  LES  equations.  This  method  was  tested  in  a 
plane  channel,  square  duct,  and  rotating  channel  by  Balaras  et  al  [8]  and  in  a  plane  channel  and 
backward- facing  step  by  Cabot  fe  Moin  [13].  More  recently,  Wang  &  Moin  [42]  used  this  method 
with  a  dynamically  adjusted  RANS  eddy  viscosity  to  perform  an  LES  of  an  airfoil  trailing  edge  flow 
[41]  with  incipient  separation.  The  results  are  better  than  those  of  the  algebraic  models,  since  the 
TBL  equations  can  account  for  more  of  the  physics  of  the  flow.  Nonetheless,  there  is  insufficient 
evidence  of  robustness  of  this  approach,  particularly  on  coarse  meshes  and  at  very  high  Reynolds 
numbers. 

To  further  assess  the  accuracy  of  this  approach  for  high  Reynolds  number  flows  under  coarse 
mesh  resolution,  a  more  severe  test  case  of  the  flow  over  a  cylinder  at  supercritical  Reynolds  number 
has  been  considered.  This  is  described  in  detail  in  Chapter  2.  The  results  are  mixed.  While  the 
method  is  capable  of  capturing  correctly  the  delayed  boundary-layer  separation  and  reduced  drag 
coefficients  after  the  drag  crisis,  the  Reynolds-number  dependence  of  the  drag  coefficient  has  not 
been  captured,  and  the  solutions  become  increasingly  inaccurate  at  higher  Reynolds  numbers.  It 
is  argued  that  at  such  marginal  resolution,  the  SGS  modeling  errors  and  numerical  errors  tend 
to  dominate  the  LES  in  the  near-wall  region,  which  cannot  be  corrected  by  a  physical  based  wall 
model. 

The  difficulty  of  formulating  a  robust  wall  model  was  highlighted  in  [12].  In  that  work,  a 
backward  facing  step  LES  was  performed  using  the  “exact”  time  series  of  the  wall  stress  from  a 
resolved  LES  as  the  wall  model.  The  results  of  this  approach  were  not  satisfactory  and  in  fact 
not  an  improvement  over  the  other  types  of  wall  models  previously  mentioned.  This  indicates  that 
SGS  and  numerical  errors  are  playing  an  important  role  in  the  coarse  grid  LES,  which  have  not 
been  accounted  for  by  the  previous  wall  models.  To  investigate  this  hypothesis  and  determine 
what  information  a  wall  model  must  provide  to  the  LES,  we  used  optimal  control  techniques  to 
compute  the  wall  stresses  in  a  channel  LES  at  Rer  —  4000.  This  work  is  discussed  in  Chapter  3 
(see  also  [33,  5]).  A  cost  function  is  defined  to  be  the  difference  between  the  plane- averaged  LES 
streamwise  and  spanwise  velocity  fields  and  their  known  mean  values  (log-law  in  the  streamwise 
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direction  and  zero  in  the  spanwise  direction).  Adjoint  equations  are  used  to  determine  the  cost 
function  derivatives,  and  iterations  are  performed  at  each  time  step  to  determine  the  best  wall 
stress.  The  simulation  results  confirm  that,  to  obtain  the  correct  mean  velocity  profile,  the  wall 
stresses  must  not  only  model  the  physics  but  also  compensate  for  numerical  and  SGS  modeling 
errors.  The  data  generated  by  this  suboptimal  control  strategy  is  then  used  to  derive  a  linear 
stochastic  estimate  (LSE)  model  which  predicts  the  wall  stresses  from  the  near-wall  velocity  field. 
When  used  in  actual  simulations,  the  LSE  model  is  shown  to  yield  accurate  mean  velocity  profiles 
over  a  range  of  Reynolds  numbers.  However,  it  is  found  to  be  sensitive  to  the  numerical  method, 
SGS  model,  and  grid  employed. 

While  valuable  insight  has  been  gained  from  the  above  control-based  wall  modeling,  the  method 
is  not  useful  as  a  predictive  tool  since  a  target  velocity  profile  must  be  specified  a  priori .  The 
computational  cost  for  the  sub-optimal  control  is  very  high  since  it  requires  both  the  solution 
of  adjoint  equations  and  many  iterations  to  achieve  convergence.  Furthermore,  the  LSE  models 
generated  from  such  computations  are  too  sensitive  to  the  numerical  parameters  to  construct  a 
universal  LSE  coefficient  database.  To  overcome  these  drawbacks,  we  have  explored  a  low-cost 
modeling  approach  using  techniques  borrowed  from  optimal  shape  design  [31]  and  a  RANS  model 
to  formulate  a  feedback  controller.  The  mathematical  framework  and  some  test  results  are  presented 
in  Chapter  4. 
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Chapter  2 

Wall  Modeling  in  High  Reynolds 
Number  Flow  over  a  Cylinder 


As  mentioned  in  the  previous  chapter,  wall  models  based  on  TBL  equations  and  their  simplified 
forms  [8,  13,  42]  have  received  much  attention  in  recent  years.  These  models,  used  with  a  RANS 
type  of  eddy  viscosity,  have  shown  promise  for  complex-flow  predictions.  The  objective  of  this 
chapter  is  to  further  assess  the  viability  and  accuracy  of  this  technique  for  high  Reynolds  number 
complex  wall-bounded  flows. 

The  flow  around  a  circular  cylinder  represents  a  canonical  problem  for  validating  new  approaches 
in  computational  fluid  dynamics.  It  is  therefore  reasonable  or  even  necessary  to  subject  the  LES 
wall  modeling  methodology  to  the  same  “grand  challenge”.  To  take  the  best  advantage  of  wall 
modeling,  we  concentrate  on  the  super-critical  flow  regime  in  which  the  boundary  layer  on  the 
cylinder  becomes  turbulent  prior  to  separation.  This  is,  to  our  knowledge,  the  first  such  attempt 
using  LES,  although  a  related  method  known  as  detached-eddy  simulation  (DES),  in  which  the 
entire  attached  boundary  layer  is  modeled,  has  been  tested  in  this  type  of  flow  [39].  Breuer  [11] 
recently  conducted  an  LES  study  at  a  high  sub-critical  Reynolds  number  of  Reo  =  1.4  x  105,  and 
showed  fairly  good  comparison  with  experimental  data  in  the  near  wake.  In  the  present  work,  three 
simulations,  at  Rep  =  5  x  105,  1  x  106,  and  2  x  106,  have  been  performed.  Simulation  results  and 
comparisons  with  experimental  data  are  summarized  below. 

2.1  Numerical  method  and  procedure 

The  same  LES  code  and  wall  model  implementation  as  used  in  [42]  are  used  for  the  present  cal¬ 
culations.  The  energy-conservative  numerical  scheme  is  of  hybrid  finite- difference/spectral  type, 
written  for  a  C-mesh  [28].  The  time  advancement  is  achieved  by  the  fractional-step  method,  in  com¬ 
bination  with  the  Crank-Nicolson  method  for  viscous  terms  and  third-order  Runge-Kutta  scheme 
for  convective  terms.  A  multi-grid  iterative  procedure  is  used  to  solve  the  Poisson  equation  for 
pressure.  The  SGS  stress  tensor  is  modeled  using  the  dynamic  SGS  model  [17,  25]. 

The  computational  domain  has  a  spanwise  size  of  2D  ( D  =  cylinder  diameter),  over  which  the 
flow  is  assumed  periodic  and  48  grid  points  are  distributed  uniformly.  In  the  planes  perpendicular 
to  the  span,  401  x  120  grid  points  are  used  in  the  C-mesh,  extending  approximately  22 D  upstream 
of  the  cylinder,  17 D  downstream  of  the  cylinder,  and  24D  into  the  far-field.  Potential-flow  solutions 
are  imposed  as  boundary  conditions  in  the  far-field,  and  convective  boundary  conditions  are  used 
at  the  outflow  boundary.  Running  at  a  maximum  CFL  number  of  1.5,  the  non-dimensional  time 
step  AtUoo/D  typically  varies  between  0.0030  and  0.0045.  To  obtain  the  results  presented  here, 
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Figure  2.1:  Instantaneous  vorticity  magnitude  at  a  given  spanwise  cut  for  flow  over  a  circular 
cylinder  at  Ren  =  106.  25  contour  levels  from  uD/Uqo  =  1  to  ljD/Uoo  =  575  (exponential 
distribution)  are  plotted. 


the  simulations  have  advanced  at  least  150  dimensionless  time  units.  The  statistics  are  collected 
over  the  last  75  or  so  time  units. 

Approximate  boundary  conditions  on  the  cylinder  surface  are  imposed  in  terms  of  wall  shear 
stress  estimated  from  a  wall  model  of  the  form 


dui 

+  vt)  -w- 

dy 


1  dp 
pdxC 


i  =  1,3 


(2.1) 


This  is  a  simpler  variant  of  the  TBL  equation  model  (cf.  Eq.  (1.2))  which  allows  for  easier  imple¬ 
mentation  and  lower  computational  cost.  Although  Wang  &  Moin  have  shown  [42]  that  the  full 
TBL  equations  (with  dynamically  adjusted  ist)  give  better  results  in  their  trailing-edge  flow,  the 
discrepancy  may  be  partly  related  to  a  surface  curvature  discontinuity  which  is  absent  from  the 
cylinder  surface.  Since  the  pressure  is  taken  from  the  LES  at  the  edge  of  the  wall  layer,  Eq.  (2.1) 
can  be  integrated  to  the  wall  to  obtain  a  closed-form  solution  for  the  wall  shear  stress  components 
[42] 


T~wi 


p  f  _I_dp  fs  ydy  1 

pdxij0v  +  ut] 


(2.2) 


where  u$i  denotes  the  tangential  velocity  components  from  LES  at  the  first  off- wall  velocity  nodes, 
at  distance  8  from  the  wall.  In  attached  flows  these  nodes  axe  generally  placed  within  the  lower 
edge  of  the  logarithmic  layer.  In  the  present  flow,  however,  8 +  (in  wall  units)  is  found  to  vary 
from  0  to  100  depending  on  the  local  skin  friction.  The  eddy  viscosity  is  modeled  by  a  damped 

mixing-length  model:  vt/v  —  ^1  —  e~y™!A^j  ,  where  =  ywurjv  is  the  distance  to  the  wall 

in  wall  units,  k  =  0.4,  and  A  =  19. 


2.2  Results  and  discussion 

In  Figure  2.1,  the  contours  of  the  vorticity  magnitude  at  a  given  time  instant  and  spanwise  plane 
are  plotted  for  Rep  =  106.  Large  coherent  structures  are  visible  in  the  wake,  but  they  are  not 
as  well  organized  and  periodic  as  in  typical  Karman  streets  at  lower  (sub-critical)  and  higher 
(post-critical)  Reynolds  numbers.  Compared  to  flows  at  lower  Reynolds  number  (e.g.  [23,  11]),  the 
boundary-layer  separation  is  much  delayed  and  the  wake  is  narrower,  resulting  in  a  much  smaller 
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Figure  2.2:  Mean  pressure  distribution  on  the  circular  cylinder. - ,  Present  LES  at  Re&  =  106; 

o  ,  Experiment  of  Warschauer  &  Leene  [43]  at  Rep  —  1.26  xlO6  (spanwise  averaged);  A  ,  Experiment 
of  Flachsbart  (in  [45])  at  Reo  =  6.7  x  105. 


drag  coefficient.  Note  that  the  rather  thick  layer  seen  along  the  cylinder  surface  consists  mostly  of 
vorticity  contours  of  small  magnitude.  These  levels  are  necessary  for  visualizing  the  wake  structure, 
but  are  not  representative  of  the  boundary-layer  thickness.  The  true  boundary  layer,  with  strong 
vorticity,  is  extremely  thin  in  the  attached  region. 

A  comparison  with  two  sets  of  experimental  data  of  the  mean  pressure  distribution  on  the 
cylinder  surface  is  depicted  in  Figure  2.2.  Very  good  agreement  is  observed  between  the  LES 
at  Re d  —  106  and  the  experiment  of  Warschauer  &  Leene  [43]  which  was  performed  at  i?e#  = 
1.26  x  106.  The  original  Cp  data  in  [43]  exhibit  some  spanwise  variations;  for  the  purpose  of 
comparison  the  average  value  is  plotted.  Relative  to  the  measurements  of  Flachsbart  (see  [45])  at 
Reo  —  6.7  x  105,  the  LES  Cp  shows  smaller  values  in  the  base  region.  Note  that  Flachsbart’s 
data  contain  a  kink  near  9  =  110°,  indicating  the  presence  of  a  separation  bubble.  This  type  of 
separation  bubble  is  characteristic  of  the  critical  regime,  and  is  difficult  to  reproduce  experimentally 
or  numerically  due  to  sensitivity  to  disturbances. 

In  Thble  2.1,  we  compare  the  mean  drag  coefficient,  the  base  pressure  coefficient,  and  the  Strouhal 
number  from  the  LES  at  Rep  =  106  with  the  experimental  values.  The  agreement  with  the  mea¬ 
surements  of  Shih  et  al  [38]  is  reasonably  good.  The  LES  somewhat  overpredicts  the  drag  coefficient 
compared  with  Shih  et  al  [38],  but  underpredicts  it  relative  to  Achenbach  [2]  (cf.  Figure  2.3).  The 
Strouhal  number  of  0.22  from  Shih  et  al  is  for  a  rough-surface  cylinder;  no  coherent  vortex  shedding 
was  observed  for  smooth  cylinders  at  Rep  larger  than  4  x  105.  Indeed,  it  is  generally  accepted  that 
periodic  vortex  shedding  does  not  exist  in  the  super-critical  regime  of  flow  over  a  smooth  cylinder 
[45].  From  our  simulation,  a  broad  spectral  peak  of  the  unsteady  lift  centered  at  St  «  0.28  is  found. 
It  can  be  argued  that  although  the  LES  is  performed  for  a  smooth  cylinder,  the  discretization  of 
the  cylinder  surface  and  the  numerical  errors  due  to  under-resolution  may  act  as  equivalent  surface 
roughness,  causing  the  flow  field  to  acquire  some  rough- cylinder  characteristics.  The  flow  at  high 
Reynolds  number  is  very  sensitive  to  surface  roughness  and  to  the  level  of  free-stream  turbulence, 
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Table  2.1:  Drag,  base  pressure  coefficient  and  Strouhal  number  for  the  flow  around  a  circular 
cylinder  at  a  Reynolds  number  of  106. _ 


Cd 

Cpfiase 

St 

"LES 

0.31 

0.32 

0.28 

Exp.  (Shih  et  al.  [38]) 

0.24 

0.33 

0.22 

Exp.  (Others,  see  [45]) 

0.17-0.40 

- 

0.18-0.50 

Figure  2.3:  Drag  coefficient  as  a  function  of  Reynolds  number.  - ,  Achenbach  [2];  •  ,  Present 

LES. 

which  contribute  to  the  wide  scatter  of  Cd  and  St  among  various  experiments  in  the  literature 
[45],  listed  at  the  bottom  of  Table  2.1.  Other  factors  causing  the  data  scatter  include  wind-tunnel 
blockage  and  end-plate  effects.  Our  simulation  results  fall  easily  within  the  experimental  range. 
Generally  speaking,  there  is  a  lack  of  detailed  experimental  data  at  super-critical  Reynolds  num¬ 
bers.  In  particular,  velocity  and  Reynolds-stress  profile  measurements  are  non-existent,  making  a 
more  detailed  comparison  impossible. 

To  assess  the  robustness  of  the  computational  method,  we  have  performed  simulations  at  Reo  = 
5  x  105  and  2  x  106,  in  addition  to  the  initial  attempt  at  Rep  =  1  x  106.  The  predicted  mean 
drag  coefficients  are  plotted  in  Figure  2.3  along  with  the  drag  curve  of  Achenbach  [2].  While  the 
simulations  predict  Cd  rather  well  at  the  two  lower  Reynolds  numbers,  the  discrepancy  becomes 
large  at  Rcd  =  2  x  106.  More  significantly,  the  LES  solutions  show  relative  insensitivity  to  the 
Reynolds  number,  in  contrast  to  the  experimental  data  which  exhibit  an  increase  in  Cd  with 
Reynolds  number  after  the  drag  crisis.  Similar  Reynolds-number  insensitivity  has  been  observed 
for  the  other  quantities  shown  previously. 

Finally,  the  skin-friction  coefficients  predicted  by  the  wall  model  in  the  LES  calculations  axe 
plotted  in  Figure  2.4  against  the  experimental  data  of  Achenbach  [2]  at  Ren  =  3.6  x  106.  The 
levels  are  very  different  on  the  front  half  of  the  cylinder,  but  are  in  reasonable  agreement  on  the  back 
half.  The  boundary-layer  separation  and  the  recirculation  region  are  captured  rather  well  by  the 
LES,  indicating  that  they  axe  not  strongly  affected  by  the  upstream  errors.  The  different  Reynolds 
numbers  in  the  LES  and  the  experiment  can  account  for  only  a  small  fraction  of  the  discrepancy. 
Note  that  our  computed  Cf  values  are  comparable  to  those  reported  by  Travin  et  al.  [39]  using  DES. 
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0.015 


Figure  2.4:  Skin  friction  distribution  on  the  cylinder  from  LES:  - ,  Rep  =  5  x  105; - , 

Reo  =  1  x  106; - ,  Reo  =  2x  106.  o  ,  Experiment  of  Achenbach  [2]  at  Reo  =  3.6  x  106. 

Travin  et  al.  attribute  the  overprediction  of  Cf  before  separation  to  the  largely-laminar  boundary 
layer  in  the  experiment,  which  has  not  been  modeled  adequately  in  either  simulation.  Poor  grid 
resolution,  which  becomes  increasingly  severe  as  the  Reynolds  number  increases,  is  another  potential 
culprit  in  the  present  work.  In  addition,  an  overprediction  of  the  skin  friction  by  the  present  wall 
model  has  also  been  observed  by  Wang  &  Moin  [42]  in  the  acceleration  region  of  the  trailing- 
edge  flow,  suggesting  that  this  simplified  model  may  have  difficulty  with  strong  favorable  pressure 
gradients. 


2.3  Concluding  remarks 

A  numerical  experiment  has  been  carried  out  to  compute  the  flow  around  a  circular  cylinder  at 
supercritical  Reynolds  numbers  using  LES.  The  simulation  is  made  possible  by  the  use  of  a  simple 
wall-layer  model,  Eq.  (2.2).  The  computational  solutions  correctly  capture  the  delayed  boundary- 
layer  separation  and  reduced  drag  coefficients  consistent  with  measurements  after  the  drag  crisis. 
In  quantitative  terms,  the  mean  pressure  distributions  and  overall  drag  coefficients  are  predicted 
reasonably  well  at  Rep  =  5  x  105  and  106.  However,  the  results  are  inaccurate  at  higher  Reynolds 
numbers,  and  the  Reynolds-number  dependence  of  the  drag  coefficient  is  not  captured.  It  must  be 
pointed  out  that  the  grid  used  near  the  cylinder  surface,  particularly  before  separation,  is  quite 
coarse  judged  by  the  need  to  resolve  the  outer  boundary-layer  scales.  This  is  in  contrast  to  the 
trailing-edge  flow  case  [42]  for  which  the  model  works  well.  At  such  marginal  resolution,  the  SGS 
modeling  errors  and  numerical  errors  tend  to  dominate  the  LES  in  the  near-wall  region,  which 
cannot  be  corrected  by  a  physical  based  wall  model  such  as  the  one  employed  here.  A  control- 
based  model,  to  be  discussed  in  the  next  chapter,  can  account  for  these  errors  and  hence  may 
provide  a  better  alternative.  In  addition,  based  on  the  experimental  observation,  the  boundary 
layer  is  largely  laminar  prior  to  separation,  which  has  not  been  modeled  adequately  by  the  current 
wall-layer  model. 
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Chapter  3 

Wall  Modeling  Using  Optimal  Control 
Techniques 


This  chapter  describes  some  recent  applications  of  control  theory  to  the  wall  modeling  problem. 
As  noted  in  Chapter  1,  there  is  evidence  that  even  the  analytically  correct  wall  stresses  would 
form  a  poor  LES  wall  model.  Chapter  2  shows  that  a  wall  model  based  purely  on  physics  does  not 
produce  accurate  solutions  for  high  Reynolds  number  flows  under  coarse  grid  resolution.  Hence,  it  is 
desirable  to  formulate  a  wall  model  with  the  capability  to  compensate  for  the  additional  numerical 
and  SGS  modeling  errors.  Since  there  is  no  known  method  for  quantitatively  determining  the  exact 
errors  a  priori  (or  creating  a  high  fidelity  model  for  them  from  a  posteriori  data),  techniques  from 
optimal  control  theory  have  been  used  to  drive  the  LES  system  towards  a  desired  state.  These 
techniques  were  chosen  since  they  seek  to  minimize  a  cost  function  whose  minimization  should 
improve  the  LES  results.  In  doing  so,  wall  stresses  that  both  compensate  for  the  missing  physics 
and  the  numerical  errors  can  be  found. 

Optimal  control  is  a  branch  of  control  theory  that  deals  with  determining  the  set  of  controls 
from  some  control  space  of  £2  functions  that  minimize  a  given  cost  function,  often  denoted  as  J . 
There  has  been  and  continues  to  be  considerable  application  of  this  technique  to  problems  in  fluid 
mechanics.  Early  applications  in  conjunction  with  turbulence  simulations  involved  minimizing  drag 
in  a  channel  [15,  9].  For  a  review  of  some  of  the  more  recent  applications  of  this  approach,  see  [10] 
and  the  references  therein.  These  successes  provided  evidence  that  the  optimal  control  technique 
could  be  successfully  applied  in  situations  with  complex  physics.  A  new  challenge  in  the  context  of 
wall  modeling  is  to  use  it  in  a  situation  where  both  the  physics  and  numerics  need  to  be  accounted 
for. 

This  approach  was  motivated  by  two  results.  The  first  was  the  previously  mentioned  work  of 
Cabot  [12]  demonstrating  that  resolved  LES  wall  stresses  did  not  perform  well  as  a  wall  model.  A 
second  result  was  that  of  Bagwell  et  ai  [7,  6]  in  their  use  of  linear  stochastic  estimation  (LSE)  to 
derive  a  wall  model.  A  linear  convolution  of  form  =  Lij  *  uj (•,  ?/,  •)  with  i  =  1,3  and  j  =  1, 2, 3 
was  used  as  the  wall  model.  The  convolution  coefficients,  were  found  using  the  LSE  technique 
of  minimizing  the  mean  square  error  of  this  model  with  DNS  training  data  at  ReT  =  180  [22]. 
This  model  worked  well  in  a  channel  flow  at  the  same  Reynolds  number.  A  channel  flow  LES  was 
then  performed  using  convolution  kernels  that  were  extrapolated  to  Rer  =  640  [6].  In  this  case, 
the  LSE  model  performed  only  marginally  better  than  the  shifted  model  of  Piomelli  et  al  [36] 
(cf.  Eq.  (1.1)).  This  result  again  demonstrates  the  possibility  that  the  effect  of  the  unmodeled 
components  (numerical  and  SGS  modeling  errors)  play  a  significant  role  in  the  wall  modeling 
problem.  It  also  shows  that  if  LSE  is  to  be  used  to  develop  wall  models,  the  training  data  must 
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come  from  an  LES  at  comparable  Reynolds  numbers.  However,  these  Reynolds  numbers  are,  in 
general,  too  high  to  perform  an  LES  without  a  wall  model.  Therefore,  a  scheme  based  on  optimal 
control  techniques  will  be  used  to  perform  such  an  LES,  thereby  generating  training  data  that  is 
at  the  correct  Reynolds  number  and  will  account  for  numerical  and  SGS  modeling  errors.  The 
methodology  and  results  are  discussed  in  the  remainder  of  this  chapter.  They  are  also  described  in 
[32,  33,  5]. 

3.1  Numerical  method 

A  second-order  accurate  finite  difference  scheme  is  used  to  discretize  the  LES  equations  on  a 
staggered  grid  system  [19].  Given  the  simple  geometry  of  a  channel  flow  considered  in  this  study, 
more  accurate  (spectral)  methods  could  have  been  used.  However,  these  highly  accurate  methods 
axe  not  flexible  enough  to  handle  industrial  applications  with  complex  geometries  (e.g.  flow  around 
an  airfoil),  where  both  low-order  numerics  for  simplicity  and  wall  modeling  for  high- Reynolds 
number  boundary  layers  are  needed.  A  staggered  grid  system  is  used  to  avoid  the  decoupled 
pressure- velocity  mode  as  well  as  the  prescription  of  a  boundary  condition  for  the  pressure.  The 
time  integration  is  a  third-order  Runge-Kutta  scheme  for  all  the  convection  and  diffusive  terms. 
The  diffusive  terms  in  the  normal  direction  to  the  wall  are  not  treated  implicitly  since  only  coarse 
grid  computations  are  considered.  Periodic  conditions  are  imposed  in  the  two  directions  x\  and 
xz  (or  x  and  z)  parallel  to  the  walls  so  that  the  Poisson  equation  can  be  solved  efficiently  using  a 
FFT-based  Poisson  solver. 

The  subgrid  scale  model  is  the  Smagorinsky  model  with  the  coefficient  determined  by  the  plane- 
averaged  dynamic  procedure  [25].  Unless  otherwise  stated,  all  quantities  are  nondimensionalized 
by  the  friction  velocity,  uT)  and  channel  half-height,  h .  The  channel  walls  are  at  y  =  ±1.  The  skin 
friction  Reynolds  number  is  then  defined  as  Rer  =  uTh/v .  When  the  mean  flow  is  converged  in 
the  statistical  sense,  the  mean  streamwise  pressure  gradient  is  equal  to  the  mean  wall  stress,  that 
is,  -d(p)/dx  =  (rw)  =  1. 

Since  ‘non-resolved’  LES  is  considered  in  this  study,  the  classical  no-slip  boundary  condition  for 
the  velocity  components  is  replaced  by  a  set  of  approximate  boundary  conditions.  More  precisely, 
the  transpiration  velocity  (in  some  cases  set  to  zero)  and  the  two  shear  stresses  and  t™2  are 
provided  to  the  momentum  equation  in  the  x\—  and  x$— direction  respectively.  The  sketch  in 
Figure  3.1  shows  the  location  of  the  variables  and  boundary  conditions  in  the  staggered  grid  system. 
The  wall  normal  direction  is  X2  (or  y)  while  Ui  (or  u,v,w)  and  p  are  the  velocity  components  and 
the  pressure.  The  shear  stresses  and  t™2  are  prescribed  either  from  one  of  the  wall  stress  models 
discussed  in  Chapter  1  or  the  optimal  control  strategy.  The  wall  transpiration  velocity  is  set  to 
zero  or  determined  using  the  control  strategy. 

3.2  Optimal  formulation 

In  the  process  of  deriving  a  new  wall  model  for  LES,  it  is  crucial  to  keep  in  mind  that: 

1.  The  objective  is  to  provide  an  approximate  boundary  treatment  able  to  handle  very  large 
Reynolds  numbers.  In  this  respect,  using  DNS  data  as  a  guide  may  not  be  the  most  judicious 
choice  since  these  data  are  only  available  for  low  to  moderate  Reynolds  numbers. 

2.  In  any  coarse  grid  LES  where  the  first  grid  point  is  within  the  logarithmic  region,  the  turbulent 
integral  length  scale  (L  =  ny)  is  less  than  half  the  grid  spacing  (Ay  =  y).  As  a  consequence 
both  sub-grid  modeling  and  numerical  errors  are  important.  Therefore  the  approximate 
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First  off-wall  cell 


t12  t12  t32  t32 


Boundary  plane  y  =  +/-  h 


Figure  3.1:  Staggered  grid  system  used  in  this  study. 


boundary  condition  should  compensate  for  these  errors  if  the  correct  mean  profile  is  to  be 
obtained  in  a  coarse  grid  computation.  In  this  case  the  best  shear  stresses  to  supply  are  not 
the  physical  ones,  as  from  a  DNS. 

It  follows  that  the  reference  data  used  to  compare  or  derive  new  wall  models  should  be  obtained 
from  a  high  Reynolds  number  LES  simulation  on  a  coarse  grid.  Of  course,  such  a  simulation  requires 
a  good  model  for  the  near- wall  region  in  the  first  place.  The  optimal  control  framework  is  applied 
here  to  conduct  such  a  simulation  without  a  priori  knowledge  of  the  necessary  wall  stress  boundary 
conditions.  The  case  of  a  channel  flow  with  constant  pressure  gradient  is  considered  ( d{p)/dx  =  1). 
The  objective  is  to  optimize  the  shear  stresses  and  r^,  and  in  some  cases  the  transpiration 
velocity  vw ,  in  order  to  minimize  a  given  cost  function.  The  mathematical  formulation  is  detailed 
in  the  following  subsections.  It  is  given  for  the  most  general  case,  in  which  wall  stresses  and 
transpiration  velocity  act  as  controls,  and  penalty  terms  for  the  deviation  from  the  mean  velocity 
profile,  deviation  from  the  rms  velocity  fluctuations,  and  control  penalty  are  included  in  the  cost 
function.  Simpler  cases  may  easily  be  derived  by  setting  certain  terms  to  zero. 


3.2.1  State  equation 


The  problem  considered  is  governed  by  the  unsteady,  incompressible,  filtered  Navier-Stokes  equa¬ 
tions  as  well  as  the  divergence-free  constraint  which  arises  from  continuity.  The  governing  equations 
read: 


dui  duiUj 
dt  +  dxj 


dp 

dxi 


-f  Sn  + 


d 

dxj 


duj 

dxj 


=  0, 


(3.1) 
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where  the  Sn  represents  the  constant  streamwise  pressure  gradient.  Note  that  no  specific  notation 
is  used  to  describe  the  spatial  filtering  associated  with  the  LES  formulation.  Each  variable  in  the 
previous  and  subsequent  equations  should  be  understood  as  a  low-pass  filtered  version  of  the  actual 
variable  (e.g.  where  the  overbar  stands  for  the  filtering  operator).  Eq.  (3.1)  is  valid  for  any 

subgrid  model  based  on  the  Boussinesq  assumption.  The  boundary  conditions,  specified  in  terms 
of  the  control  parameter  (j>  defined  below,  on  Eq.  (3.1)  axe: 

du  {  dvn  _  _1_ 

dyn  dx  uw 

Vn  —  (3*2) 

dw  t  dvn  _  J_ 

dyn  dz  vjw 

where  the  subscript  n  stands  for  the  outward  normal  to  the  wall  and  vw  is  the  wall  value  of  the 
total  dynamic  viscosity  v  +  v%  (in  this  work  vw  =  u).  The  control  parameter  (f)  is  defined  as 
0=  (4>1 i,<t>v,<t>w)  =  {t?2 iVw,t£2)  at  y  =  +1  and  (j>  =  {(fru^vAw)  =  -(r$,  vw,rg>)  at  y  =  -1. 

In  the  classical  optimal  control  procedure  the  objective  is  to  reduce  the  given  cost  function  for 
some  period  of  time.  This  method  has  been  proven  to  be  efficient  [1].  However,  this  is  a  very 
expensive  procedure  in  terms  of  storage  and  manipulation  of  many  3D  fields  over  the  entire  period 
under  consideration  [15].  We  therefore  make  use  of  a  more  affordable  sub-optimal  procedure  in 
which  the  state  equation  is  first  discretized  in  time,  then  a  control  procedure  is  used  to  minimize 
the  cost  function  over  a  short  period  of  time  (the  time  step)  at  each  time  step  [9].  This  method  does 
not  necessarily  provide  the  ‘best’  answer  but  it  is  much  more  cost  effective  than  the  optimal  strategy. 
The  equation  of  state  (3.1)  is  therefore  discretized  in  time  by  assuming  an  implicit  discretization: 


dui 

dxj 


,  duA 

w 

n-f  1 

dxi) 

)\ 

—2f3At 

»»r' 

dxj 

RHS71 

(3.3) 

0 


The  boundary  conditions,  Eq.  (3.2),  apply  to  Eq.  (3.3)  and  the  second  order  Crank-Nicolson 
method  is  recovered  if  the  parameter  /?  is  chosen  to  be  /?  =  Note  that  in  Eq.  (3.3)  only  the 
terms  involving  the  solution  at  the  current  time  step  n  +  1  are  written  explicitly.  This  is  because 
in  the  sub-optimal  framework  developed  in  Bewley  et  al  [9],  only  the  terms  at  time  n  +  1  in  the 
semi-discretized  state  equation  are  used  in  the  optimization  process.  The  terms  which  depend 
only  on  the  variables  at  the  previous  time  step  n  are  gathered  in  the  generic  notation  RHSn  and 
disappear  in  the  analytical  development.  Note  that  there  is  a  slight  inconsistency  between  the  time 
integration  used  to  solve  the  state  equation  (explicit  third-order  Runge-Kutta,  Section  3.1)  and 
time  integration  supposed  in  the  process  of  defining  the  control  problem  (implicit  Crank-Nicolson). 
The  optimal  control  problem  and  the  resulting  optimization  algorithm  proposed  in  Section  3.2.4 
may  therefore  be  considered  as  only  approximation  of  the  ‘exact*  control  problem  that  should  be 
defined  and  solved  for  the  flow  problem  considered.  However,  in  view  of  the  results  shown  in 
Sections  3.3  and  3.5,  and  the  small  time  steps  used  for  the  explicit  scheme,  this  approximation  is 
sufficient  to  provide  valuable  insight  for  identifying  a  practical  wall  model  for  coarse  grid  LES. 


3.2.2  Objective  function 

In  the  sub-optimal  control  approach,  the  boundary  conditions  (specified  by  the  control  parameter 
(j>)  are  used  as  control  to  minimize  an  objective  function  at  each  time  step.  The  goal  is  to  provide 
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numerical  boundary  conditions  to  the  flow  solver  so  that  the  overall  solution  is  consistent  with  what 
is  expected  in  a  channel  flow.  The  objective  function  is  specified  as  follows: 

3  3 

J(u;  (f>)  -  ^2  ^mean,i(w;  (f>)  +  i7rms,i(^;  <£)  +  ^2  ^penalty,* WO-  (3*4) 

i=l,3  i=l  1 

The  objective  function  consists  of  the  three  components.  Jinean  measures  the  distance  from  the 
plane-averaged  LES  solution  to  a  desired  reference  velocity  profile.  The  second  component,  Jrms 
measures  the  distance  from  the  plane-averaged  velocity  fluctuation  intensities  to  desired  target  pro¬ 
files.  Finally,  the  third  component,  Jpena]ty  penalizes  the  use  of  large  controls  <j).  The  component 
objective  functions  are  defined  below. 

For  the  mean  streamwise  velocity  the  target  or  reference  profile  is  taken  as  a  logarithmic  velocity 
profile  throughout  the  channel:  —  «-1  lny+  +  C.  The  spanwise  velocity  reference  profile  is 

simply  u3  ref  =  0.  The  difference  between  the  reference  velocity  profile  and  the  plane-averaged  LES 
solution  is  a  function  of  the  wall-normal  coordinate,  y,  and  can  be  expressed  as 

K  (y)  =  jf  j K  -  Vef) dxdz  (*  =  3)  (3-5) 

where  A  is  the  channel  area  in  the  homogeneous  plane.  Note  that  any  reference  profile  suitable  for 
a  parallel  flow  could  have  been  used.  Notably,  a  more  realistic  shape  could  have  been  used  near 
the  channel  center.  However,  the  logarithmic  profile  is  well  suited  for  the  near-wall  region  since 
we  are  using  a  coarse  mesh  and  the  Reynolds  number,  Rer  =  4000,  is  sufficiently  high  so  that  the 
first  grid  point  lies  in  the  logarithmic  region  (y+  «  121).  The  mean  component  of  the  objective 
function  is  then: 

i7mean,i(«;  4>)  =  ai  J  i  Sm(y)2  dV,  (*  =  1,3)  (3.6) 

Note  that  there  is  no  need  to  specify  a  target  profile  for  the  plane-averaged  wall-normal  velocity 
since  that  will  be  identically  zero  at  each  time  step  provided  there  is  no  net  transpiration  velocity 
through  the  boundaries. 

The  velocity  fluctuation  intensities  are  targeted  through  the  Jrms  component  of  the  objective 
function.  The  plane-averaged,  mean  square  velocity  fluctuations  are  compared  at  each  time  step 
to  the  mean  square  velocity  fluctuations,  (^ref)2?  from  the  LES  of  Kravchenko  et  al.  [24]  which 
was  performed  at  the  same  Reynolds  number  using  a  zonally  defined  mesh  to  resolve  the  near¬ 
wall  region.  The  distance  between  the  plane-averaged  mean  square  velocity  fluctuations  and  their 
reference  profiles  can  be  measured  as 

*«{(y)  =  ~iff  ((Ui  ~  (ui))2  -  ui%f)  dxdz>  (*  =  3)>  (3-7) 

where  (ui)  denotes  the  average  over  the  homogeneous  directions  of  the  velocity  component  The 
velocity  fluctuation  intensity  component  of  the  objective  function  is  then 

/+i 

t  &u\{y)2dy  (*  =  1,2,3)  (3.8) 

Finally,  to  prevent  numerical  instabilities  it  is  necessary  to  regularize  the  control,  that  is,  the 
approximate  boundary  conditions,  by  including  a  penalty  component  in  the  overall  objective  func¬ 
tion: 

^penalty  ,*W)  ~  J  ±1  ^ Ui  ^X(^z  ^  J  dxdz  (i  =  1,2,3).  (3.9) 
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The  first  term  in  the  penalty  component  attempts  to  prevent  the  mean  square  norm  of  the  control 
parameter  from  becoming  too  large.  In  the  case  of  transpiration  velocity  control,  it  was  found 
necessary  to  prevent  the  transpiration  velocity  from  becoming  too  large  at  any  single  point,  hence 
the  second  term  in  the  penalty  component  (3.9)  was  added. 

Note  that  each  component  of  the  objective  function  includes  a  scalar  parameter:  a*,  /%,  7*,  or  A. 
These  scalars  allow  the  relative  importance  of  the  various  objectives  to  be  changed  in  the  overall 
objective. 


3.2.3  Adjoint  problem 

The  gradient  of  the  objective  function  J  with  respect  to  the  control  parameter  is  estimated  by 
using  the  Frechet  differential  [40]  defined  for  any  functional  F  as: 


DF~=  Fil+e^-m_ 

D<t> V  e->o  e  ’ 

where  <j)  is  an  arbitrary  direction.  Prom  Eq.  (3.4)  the  gradient  J  is: 
DJ- 


D(j) 


<f>  =  Y''  oti  f  [  [  2^j-Uidxdydz 

*=i,3  J  '  Jn 


Ui  dx  dy  dz 


(3.10) 


(3.11) 


+A-^  J  J  ^4>l2j>U2dxdz  +  J  J  ^KiKidxdz 


where  Ui  denotes  the  Frechet  derivative  of  Uj1 .  The  gradient  of  J  cannot  be  calculated  directly 
from  Eq.  (3.11)  since  the  derivatives  Ui  are  unknown. 

We  now  formulate  an  adjoint  problem  to  find  the  gradient  (3.11).  The  first  step  is  to  take  the 
derivative  of  Eq.  (3.3)  with  respect  to  the  control  (j) : 


with  boundary  conditions: 


dU_  dVn  ^7 

dyn  dx  vju 

Vn  ~  (fiv 

dw  d}^  _ 

dyn  dz  vw  w 


0 

0 


(3.12) 


(3.13) 


The  right-hand  side  term  in  Eq.  (3.12)  is  now  zero  since  the  flow  field  at  time  step  n  does  not 
depend  on  the  control  4>  for  the  current  time  step.  Therefore,  the  superscript  ‘n  +  1’  has  been 

S  t  (it|  —  (in)) 

Technically,  the  second  term  in  Eq.  (3.11)  should  include  the  integral  f  f  Jn  2-^ — ^ - (Ui  —  (Ui))  dx  dy  dz,  but 

we  make  the  approximation  that  (Ui)  =0  since  | (£4 ) |  <  \Ui\  in  general. 
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dropped  for  clarity.  Note  also  in  Eq.  (3.12)  that  the  Frechet  derivative  of  the  eddy  viscosity  was 
assumed  to  be  zero,  that  is,  Dvt/D(f)  —  0.  The  latter  approximation  can  be  justified  for  short  time 
intervals  [14].  Moreover,  this  system  of  equations  is  linear  in  the  variables  Ui  and  P,  where  V  is 
the  Frechet  derivative  of  the  pressure.  Therefore  it  can  be  written  in  the  form: 

AQ  =  0,  (3.14) 

where  A  is  the  linear  operator  acting  on  the  vector  0  =  (UiJV)'r.  The  linear  system  (3.14)  with 
unknown  boundary  conditions  (3.13)  cannot  be  solved  directly;  instead,  an  adjoint  operator,  A *, 
is  formulated  by  considering  the  equation 

(AG,*)  =  (0,^4*^)  -f  BT,  (3.15) 

where  (•,•)  stands  for  the  inner  product  defined  as  the  integral  over  the  flow  domain  of  the  dot 
product  of  the  two  vectors  and  4/  is  the  adjoint  state  vector  =  (rji,7v)T.  Finding  the  adjoint 
operator,  A *,  and  the  boundary  terms,  BT,  is  a  straightforward  exercise  in  integrating  by  parts. 
The  adjoint  operator  acting  on  the  adjoint  state  vector,  that  is,  is  defined  by  the  equations: 


f  ^  +  At  itt+  ~  ui^x 

~  ~  ^7  (("  +  vt)  +  ^))_ 

i 

> 

i§ 

A*V 

and  the  boundary  terms  are: 

BT=A,/t 


±i 


(Press  +  Conv  4-  Vise)  dx  dz 


with 


Press 

Conv 


Vrj2n  -  Vn7r 
7] iLfiVn 


-  -  -[*(£*£)-“(£*» 

Prom  Eq.  (3.14),  the  relation  (3.15)  defining  the  adjoint  operator  reduces  to 

(A*V,e)  =  -BT. 

3.2.4  Gradient  estimate 


(3.16) 


(3.17) 


(3.18) 


(3.19) 


We  now  have  the  liberty  to  choose  boundary  conditions  and  right-hand  side  terms  for  the  adjoint 
problem  such  that  the  relation  (3.19)  can  be  utilized  to  calculate  the  gradient  of  J .  By  comparing 
Eqs.  (3.11),  (3.17),  (3.18)  and  (3.19),  it  appears  that  a  judicious  choice  for  the  definition  of  the 
adjoint  problem  is: 

/  ai  6U  +  fa 5U>  (u  -  (u))  \ 


A*V  = 


with  boundary  conditions  at  the  wall: 


p2Sv'{v  -  ( V )) 
a3^w  +  p3&w'(™  -  (w)) 

0 


(3.20) 


/ 


rnvn  +  vw— 
dyn 

V2n 

dm 


mVn  +  Vw 


dyn 


0 

0 

0. 


(3.21) 
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In  doing  so,  Eq.  (3.19)  can  be  re-written  as: 


= 

D<f> * 


T)l<i>u  +  I  v  —  2l/, 


-)4 

Vn  ) 


+  m<f>w  \  dxdz 


+Xj  j  J  ^l^U2dxdz  +  Yy)i^  J  f_±i^  dxdz 


Since  Eq.  (3.22)  is  valid  for  all  directions  <f>,  the  gradient  of  J  may  be  extracted: 


271 

=  AtT)l>w  + 


=  Afl  7T, 


,  272  ,  4A  3 

+  T*2  +  T*! 


(3.22) 


27^ 

=  Atrj3jW  +  -^-<£35 


where  the  subscript  w  stands  for  the  values  at  the  wall.  A  control  procedure  using  a  simple  steepest 
descent  algorithm  at  each  time  step  may  now  be  proposed  such  that: 


±n+l,k+l  _  ±n+l,k  _  umr +1’fc) 

<p  -<p  V  D(j) 


(3.23) 


where  the  parameter  \i  can  be  varied  to  change  the  rate  of  convergence  and  the  extra  superscript 
k  refers  to  the  subiterations  in  the  descent  algorithm.  Note  that  the  adjoint  operator  depends  on 
the  state  vector  {ui,P)T  at  time  n  -b  1  so  that  the  state  equation  and  the  adjoint  problem  must  be 
solved  simultaneously  to  obtain  the  sub-optimal  approximate  boundary  conditions.  The  following 
algorithm  is  used  a  each  time  step  of  the  flow  solver  to  obtain  the  optimized  boundary  conditions: 

1.  Start  with  the  state  vector  (u*,  P)T  at  time  n,  the  adjoint  vector  (77*,  7r)r  and  control  parameter 
cf)  at  sub-iteration  n*, 

2.  Use  Eq.  (3.20)  with  boundary  condition  (3.21)  to  compute  the  adjoint  vector  at  sub-iteration 
nfc+i.  For  this  purpose,  the  operator  for  the  adjoint  velocity,  Eq.  (3.16),  is  re-written  as: 


-77nfc+1  +  - 
2k  2 


h + **  [£  - +  - ~  w,  (<— )  (§  - +  fe) )]  <-> 


where  only  the  first  term  is  taken  at  sub- iteration  n*+i,  the  others  being  considered  at  sub¬ 
iteration  n^.  A  Poisson  equation  is  solved  for  the  adjoint  pressure  at  each  sub-iteration  to 
enforce  the  divergence-free  condition  on  77”** 1 . 

3.  Compute  the  gradient  of  the  cost  function  at  sub-iteration  1  by  using  Eq.  (3.22), 

4.  Compute  the  control  parameter  4>  at  sub- iteration  by  using  Eq.  (3.23), 

5.  Compute  the  £2-norm  of  the  difference  0n*+l  —  </>nk. 


6.  If  it  is  more  than  0.1%  of  the  jC2-norm  of  then  go  back  to  step  1.  Otherwise,  use  (f)u  and 
4>w  as  approximate  boundary  conditions  to  compute  the  state  vector  at  time  n  +  1. 


16 


3.2.5  Validating  the  gradient  of  the  objective  function 

To  validate  the  gradient  computation,  finite  difference  approximations  to  the  gradient  were  calcu¬ 
lated.  This  is  relatively  simple  to  do.  Given  a  control  vector  0  and  a  velocity  field  w,  choose  a 
small  value  of  e  and  perturb  the  control  vector  at  one  point  by  the  amount  e  (e.g.  add  e  to  at 
one  point  on  the  lower  wall)  to  obtain  a  new  control  vector  0  +  e0.  Now  advance  the  velocity  field 
one  time  step  and  explicitly  calculate  the  value  of  the  objective  function  (3.4),  that  is,  calculate 
e0).  The  approximate  gradient  in  the  direction  0  is  then: 

-=— 0  ~  - 

D(j)  € 

By  comparing  the  approximation  (3.25)  to  a  centered  difference  approximation,  it  was  found 
that  e  =  10“3  produces  good  approximations  to  the  gradient.  By  successively  perturbing  the 
control  vector  0  at  every  point,  it  is  possible  to  approximate  the  entire  gradient  DJ/Dcf).  This 
finite  difference  gradient  approximation  can  then  be  compared  to  the  gradient  approximated  by 
the  adjoint  problem  described  above. 

It  was  found  that  the  correlation  between  the  two  gradient  approximations  was  generally  in 
excess  of  80%.  When  the  transpiration  velocity  was  not  included  in  the  control,  the  correlation  was 
generally  in  excess  of  90%.  Thus  we  are  led  to  believe  that  the  adjoint  problem  defined  above  may 
be  yielding  satisfactory  approximations  to  the  gradient  of  the  objective  function,  but  further  work 
is  necessary  to  determine  if  the  gradient  approximation  can  be  improved. 

3.3  Results:  optimal  control  without  transpiration  velocity 

Several  LES’s  have  been  performed  to  test  the  optimal  control  strategy  described  in  the  previous 
section.  We  first  investigate  the  case  when  only  the  wall  stress  boundary  conditions  are  used  as 
control.  The  wall  normal  velocity  is  set  to  zero  at  the  boundary:  vn  =  0V  =  0.  In  addition,  we 
consider  an  objective  function  for  mean  flow  only,  by  setting  /%  to  zero.  The  a^s  are  taken  such 
that  a\  =  <23  =  1.  Several  numerical  tests  were  performed  to  fix  the  coefficients  //,  71,  and  73.  It 
was  found  that  the  value  /z  =  5  x  103  ensured  good  convergence  of  the  steepest  descent  algorithm, 
while  71  =  73  =  4  x  10~5  ensured  that  the  whole  algorithm  is  stable. 

3.3.1  Statistics 

In  this  section,  the  grid  is  uniform  in  all  directions  with  32  x  33  x  32  cells  and  the  domain  size  is 
(27rh,  2/i,  2irh/3)  where  h  is  the  channel  half  height.  The  Reynolds  number  based  on  the  friction 
velocity  ur  and  h  is  4000.  In  wall  units,  the  grid  spacing  is  Ax+  «  785,  At/+  «  242  and  A z+  «  262. 
Since  a  staggered  mesh  system  is  used,  the  first  u  velocity  point  is  located  at  y+  «  121.  At  this 
coarse  resolution,  most  of  near-wall  turbulent  structures  cannot  be  captured  so  that  an  effective 
wall  model  is  necessary  to  compensate.  The  computation  was  run  for  a  sufficiently  long  period 
to  be  statistically  independent  of  the  initial  condition  and  then  statistics  were  accumulated  over  a 
time  period  of  order  20  h/uT. 

Figure  3.2  shows  the  mean  velocity  profile  from  the  LES  in  which  the  optimal  procedure  of 
Section  3.2  was  used  to  obtain  the  approximate  boundary  condition.  The  mean  value  of  r\ \  was 
either  provided  by  the  optimal  procedure  itself  or  re-computed  so  that  the  first  point  coincides  with 
the  logarithmic  law  (u)+  =  2.41  In  +  5.2.  The  results  are  very  similar  in  the  two  cases.  The 
overall  agreement  is  much  better  than  with  the  shifted  model  [36],  Eq.  (1.1),  which  is  also  shown 
in  Figure  3.2  for  easier  comparison.  An  artificial  boundary  layer  still  develops  between  the  second 


(3.25) 
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Figure  3.2:  Mean  velocity  profiles  for  Rer  =  4000.  -  ,  control  includes  wall  stresses  and 

transpiration; - ,  control  includes  wall  stresses  only; - ,  no  control,  uses  wall  stress  model 

of  [36]; .  ,  logarithmic  reference  profile,  =  2.41  In  j/+  +  5.2. 


and  the  third  grid  point  but  its  amplitude  is  much  smaller  than  with  the  shifted  model.  The  deficit 
in  the  log- law  intercept  is  of  order  0.25  compared  to  1.35  in  [36]. 

Note  that  the  optimized  wall  stress  boundary  conditions  produce  a  mean  velocity  profile  that 
is  nearly  exact  for  the  first  two  grid  points.  The  small  error  in  the  channel  interior  is  believed 
to  be  due  to  the  sub-optimal  formulation  in  which  the  wall  stresses  axe  optimized  only  over  each 
time  step.  However,  the  result  in  Figure  3.2  shows  a  clear  improvement  in  comparison  with  the 
shifted  model.  An  additional  computation  was  performed  where  only  rf2  was  optimized  (^3  was  set 
to  zero  at  the  wall),  resulting  in  a  mean  velocity  profile  is  intermediate  between  the  two  previous 
computations.  The  conclusion  is  that  both  and  r™2  must  be  optimized  (or  modeled  carefully) 
if  a  velocity  profile  close  to  the  target  one  is  sought.  The  case  of  control  with  wall  stresses  and 
transpiration  velocity  is  also  plotted  in  Figure  3.2,  which  will  be  discussed  later  in  Section  3.4. 

3.3.2  Shear  stress  structure 

From  Eq.  (1.1),  the  wall  shear  stresses  in  the  LES  with  the  shifted  model  should  be  perfectly 
correlated  with  the  velocity  components  at  the  first  plane,  shifted  in  the  upstream  direction  by  the 
amount  As.  This  is  shown  in  Figure  3.3  which  displays  typical  iso-lines  of  and  w.  The 

upstream  shift  is  only  a  fraction  of  the  cell  size,  viz.  A5  «  0.67Ax,  and  is  hardly  visible  in  the  figure. 
Both  t\ 2  and  u  are  characterized  by  elongated  structures  in  the  streamwise  direction  while  rf2  and 
w  are  characterized  by  more  rounded  structures  in  the  x—z  plane.  The  reference  data  from  the  sub- 
optimal  computation  are  plotted  in  Figure  3.4.  The  perfect  correlation  between  the  wall  stresses 
and  the  velocity  components  at  the  first  off- wall  plane  no  longer  holds.  The  streamwise  velocity 
u  still  displays  elongated  structures  in  x  as  well  as  rf2.  In  contrast,  and  w  are  less  elongated. 
Eight  fields  from  the  sub-optimal  LES  have  been  used  to  compute  correlation  coefficients  between 
the  velocity  gradients,  the  velocity  components  at  the  first  off-wall  plane  and  the  wall  stresses. 
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Figure  3.3:  From  top  to  bottom:  Contours  of  u,  t™2  and  w  from  the  LES  with  the  shifted  model 
of  Piomelli  et  al  [36]  at  Rer  =  4000.  The  velocity  components  are  from  the  first  off- wall  plane. 


The  best  results  (those  with  the  correlation  greater  than  0.3)  are  reported  in  Figure  3.5.  The 
position  of  the  perfect  correlations  between  and  u,  and  w  which  are  assumed  in  the  shifted 
model  [36],  Eq.  (1.1),  is  also  reported.  Note  that  in  the  sub-optimal  calculation,  is  almost  not 
correlated  with  the  streamwise  velocity  downstream  shifted  by  As.  The  correlation  is  better  when 
an  upstream  shifted  version  of  u  is  used  instead.  On  the  other  hand  the  maximum  of  correlation 
between  rf2  and  w  is  located  close  to  the  assumed  downstream  shift  As.  The  best  correlation  for 
t\ 2  is  found  with  the  spanwise  derivative  of  w ,  whereas  r™2  best  correlates  with  du/dz.  Other  non 
negligible  correlation  coefficients  are  found  between  rf2  and  du/dx ,  and  r™2  and  dv/dz .  No  clear 
picture  is  available  yet  to  explain  these  correlations.  Finally,  note  the  fairly  good  and  somewhat 
expected  negative  correlation  between  and  v.  This  correlation  supports  the  underlying  idea 
in  the  ejection  model  of  Piomelli  et  al  [36].  A  more  systematic  way  of  exploiting  correlations 
between  the  approximate  boundary  condition  and  the  velocity  field  close  to  the  wall  is  presented 
in  Section  3.5. 

3.3.3  Discussion 

The  reduced  deficit  in  the  intercept  with  the  sub-optimal  computation  is  associated  with  a  better 
representation  of  the  gradient  of  the  mean  velocity  profile  within  the  first  few  grid  points.  This  is 
better  seen  in  Figure  3.6  which  displays  the  quantity  ny+  d  (u)+  /dy+  as  a  function  of  the  distance 
to  the  wall.  Theoretically,  this  quantity  should  be  equal  to  unity  in  the  case  where  the  sub- 
optimal  strategy  is  used  since  the  value  of  k  is  taken  to  be  k  =  1/2.41,  consistently  with  the  target 
velocity  profile  (u)+  =  2.41  In y+  4*  5.2.  For  the  case  where  the  shifted  model  [36],  Eq.  (1.1),  is 
used,  this  quantity  is  not  expected  to  be  exactly  unity  since  the  value  picked  for  k  is  somewhat 
arbitrary.  However  this  quantity  is  constant  if  the  mean  profile  follows  a  logarithmic  behavior. 
Practically,  ny*  d(u)+  fdy+  is  not  found  to  be  equal  to  unity  when  computed  from  the  exact  log 
profile  (u)+  =  2.41  In y+  4-5.2  because  of  the  large  errors  associated  with  the  second-order  finite 
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Figure  3.4:  From  top  to  bottom:  Contours  of  u,  r™2  and  w  from  the  LES  with  the  sub-optimal 
strategy  at  Rer  —  4000.  The  velocity  components  are  from  the  first  off- wall  plane. 
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Figure  3.5:  Best  correlation  coefficients  for  (a)  r\ 2,  (b)  r^.  - ,  sub-optimal  computation;  •  , 

shows  the  position  where  the  correlation  between  and  u,  and  rf2  and  w  in  the  shifted  model  of 
Piomelli  et  al.  [36]  is  assumed  to  be  unity. 
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Figure  3.6:  Non-dimensionalized  mean  velocity  gradient  d  (u)+  /dy+  for  Rer  =  4000.  - , 

sub-optimal  computation;  — i —  ,  Shifted  model  of  Piomelli  et  al.  [36]; - ,  not  exact  differenti¬ 

ation  applied  to  (u)+  =  2.41  In  y+  +  5.2. 


differences  on  the  coarse  grid  considered.  This  is  shown  in  Figure  3.6  which  also  demonstrates  that 
the  non-dimensionalized  mean  gradient  from  the  sub-optimal  computation  follows  closely  (to  within 
8%)  its  value  from  the  exact  log  law.  The  mean  gradient  from  the  computation  with  the  shifted 
model  [36]  is  25%  below  its  expected  value,  which  explains  the  deficit  in  the  intercept  observed  in 
Figure  3.2. 

The  mean  velocity  gradient  appears  in  the  mean  momentum  equation  which  reads: 


/  /  n  dlu) 

(uv)  =  y  +  + 

dy 


(3.26) 


This  equation  reduces  to: 


d(«)  w  (mV)  -  y 

dy  ~  v  +  (vt) 


(3.27) 


under  the  usually  well  verified  assumption  (vt^J  ~  (vt)  Eq.  (3.27)  reveals  that  for  a  given 

amount  of  eddy  viscosity,  the  mean  velocity  gradient  is  directly  related  to  the  difference  between 
the  total  (— y)  and  resolved  (—  (u'v1))  shear  stress,  viz.  (u'v')  —  y,  the  larger  the  difference  the  larger 
the  gradient.  This  is  confirmed  in  Figure  3.7  which  shows  that  the  artificial  condition  provided  by 
the  sub-optimal  strategy  leads  to  an  equilibrium  in  which  the  resolved  part  of  the  stress  is  smaller. 
It  is  also  worth  noting  that  the  quantity  (u'v')  —  y  is  quite  different  in  the  two  computations 
but  that  the  absolute  value  of  the  resolved  stress  are  both  close  to  unity  (|  (u'v1)  |max  »  0.9  for 
the  computation  with  the  shifted  model  of  Piomelli  et  al.  [36],  |  (u'v')  |max  «  0.86  for  the  sub- 
optimal  case).  Since  there  is  a  large  difference  in  the  velocity  gradient  (it  is  30%  higher  in  the 
sub-optimal  case),  the  production  of  kinetic  energy  —  (u'v')  d  (u)  /dy  is  necessarily  greater  in  the 
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Figure  3.7:  Difference  between  the  total  and  the  resolved  stress  for  Rer  =  4000.  - ,  sub-optimal 

computation;  — i —  ,  Shifted  model  of  Piomelli  et  al  [36]. 


sub-optimal  case,  as  shown  in  Figure  3.8a.  Figure  3.8b  shows  that  this  is  also  true  for  the  dissipation 
rate  of  kinetic  energy,  a  direct  consequence  of  the  fact  that  production  and  dissipation  balance 
reasonably  well  in  this  high- Reynolds  number  channel  flow.  Note  that  the  total  dissipation  is  very 
well  approximated  by  the  subgrid  scale  (SGS)  dissipation  esgs  =  (2 vtSijSij)  that  has  been  reported 
in  Figure  3.8b.  The  SGS  dissipation  can  further  be  approximated  (to  within  a  few  percents)  by 
^approx  —  2  (ut)  ( SijSij ).  The  form  of  the  approximate  dissipation  eapprox  indicates  that  the  sub¬ 
optimized  boundary  condition  can  act  on  the  flow  field  and  generate  an  equilibrium  with  higher 
dissipation  in  the  near  wall  region  by  either  increasing  the  mean  eddy- viscosity  or  increasing  the 
velocity  fluctuations.  The  mean  eddy-viscosity  profiles  in  the  two  computations  are  found  to  be 
almost  identical  (not  shown).  The  turbulent  fluctuations  are  higher  in  the  sub-optimal  case  than  in 
the  shifted  model  [36]  computation,  see  Figure  3.9.  The  agreement  with  the  fully-resolved  LES  data 
of  Kravchenko  et  al  [24]  turns  out  to  be  better  with  the  sub-optimized  boundary  condition  regarding 
the  normal  and  spanwise  direction,  worse  regarding  the  streamwise  direction.  In  any  case,  the 
reasoning  given  above  indicates  that  the  increase  in  the  turbulent  fluctuations  when  the  sub-optimal 
strategy  is  used,  Figure  3.9,  is  consistent  with  the  previous  findings  on  the  mean  velocity  gradient, 
Figure  3.2,  resolved  shear  stress,  Figure  3.7,  turbulence  production  and  dissipation,  Figure  3.8. 

Similar  results  have  been  obtained  for  two  other  Reynolds  numbers,  namely  ReT  =  640  and 
20000  (not  shown).  For  the  first  case,  the  mean  velocity  profile  from  the  sub-optimal  computation 
is  slightly  closer  to  the  logarithmic  law  than  that  obtained  from  the  shifted  model  [36].  For  the 
case  Rer  =  20000,  the  deficit  in  the  intercept  is  the  same  than  for  the  case  Rer  =  4000,  viz.  0.25. 
Also  the  profiles  of  turbulent  fluctuations  do  not  change  if  plotted  in  outer  coordinates  y/h.  This 
means  that  the  asymptotic  regime  in  terms  of  Reynolds  number  has  been  reached:  results  of  the 
same  quality  as  above  can  be  obtained  for  arbitrarily  large  Reynolds  numbers  with  the  same  coaxse 
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Figure  3.8:  Turbulence  (a)  production  (b)  subgrid  scale  dissipation  from  the  sub-optimal  compu¬ 
tation  scaled  by  the  same  quantity  from  the  LES  with  the  shifted  model  of  Piomelli  et  al.  [36]  for 
Rer  —  4000. 


Figure  3.9:  Root-mean-square  of  velocity  fluctuations  with  objective  to  control  mean  flow  only. 

,  control  includes  wall  stresses  and  transpiration; - ,  control  includes  wall  stresses  only; 

,  reference  profiles  from  [24]. 


grid. 


3.4  Results:  optimal  control  with  transpiration  velocity 

The  results  presented  in  the  previous  section  were  computed  assuming  that  the  wall- normal  velocity 
is  identically  zero  at  the  boundary.  However,  since  the  no-slip  boundary  condition  cannot  be  applied 
without  adequate  near-wall  resolution,  perhaps  it  does  not  make  sense  to  insist  that  the  velocity 
normal  to  the  boundary  is  zero.  After  all,  a  wall  model  should  capture  the  effects  of  the  near- wall 
turbulence  on  the  outer  flow,  including  such  hallmarks  of  near-wall  turbulence  as  ejections  and 
sweeps.  The  combination  of  non-zero  boundary-normal  velocity  with  wall  stresses  should  allow  the 
wall  model  to  influence  more  of  the  computational  domain  than  wall  stress  boundary  conditions 
alone,  since  the  boundary- normal  velocity  affects  the  entire  flow  directly  via  the  continuity  equation. 

3.4.1  Objective  function  for  mean  flow  only 

Transpiration  velocity  is  now  added  to  the  control  set  to  determine  if  there  is  an  improvement  of 
the  prediction  of  the  mean  velocity  profile.  To  test  the  influence  of  transpiration  velocity  only  on 
the  mean  velocity  profile,  the  constants  pi  in  the  objective  function,  Eq.  (3.4),  are  set  to  zero  so 
that  only  the  desired  mean  velocity  profile  is  targeted.  For  this  simulation  the  parameters  in  the 
objective  function  (3.4)  were:  a\  =  a3  =  l,/?i  =  P2  =  /?3  =  0,71  =  73  =  10“4,72  =  0.02,  and 
A  =  0.  The  relaxation  parameter  in  the  steepest  descent  algorithm  was  /i  =  103. 

The  new  mean  velocity  profile  is  plotted  in  Figure  3.2  as  the  solid  line,  which  shows  that,  indeed, 
the  addition  of  the  transpiration  control  improves  the  mean  profile  slightly  over  the  case  when  only 
wall  stress  controls  are  considered.  Also  shown  in  Figure  3.2  is  the  mean  velocity  profile  obtained 
by  using  the  simple  wall  stress  model  of  [36],  Eq.  (1.1),  that  correlates  the  streamwise  wall  stress  to 
the  streamwise  velocity  at  a  point  away  from  the  wall  and  slightly  downstream.  The  latter  model 
yields  results  that  are  typical  of  most  current  wall  stress  models  for  this  flow  configuration. 

The  improvement  in  the  mean  velocity  profile  is  encouraging.  However,  this  is  obtained  at  the 
expense  of  turbulence  intensity.  Figure  3.9  shows  the  root  mean  square  (rms)  velocity  fluctuations 
for  the  sub-optimal  wall  stress  boundary  conditions  with  and  without  the  addition  of  transpiration 
velocity  control.  The  rms  velocity  fluctuations  actually  increase  with  the  addition  of  transpiration, 
which  is  certainly  in  the  wrong  direction  since  the  fluctuation  intensities  are  already  over-predicted. 

3.4.2  Objective  function  including  mean  flow  and  rms  velocities 

In  this  case,  the  objectives  of  matching  the  rms  velocities  of  Kravchenko  et  at  [24]  are  added 
into  the  cost  function.  For  this  simulation  the  parameters  in  the  objective  function  (3.4)  were: 
ai  =  a3  =  l,Pi  =  p2  =  h  =  3  x  10‘4,7i  =  5  x  10“5,t2  =  HT3,73  =  4x  10"6  and  A  =  5  x  10"3. 
The  relaxation  parameter  in  the  steepest  descent  algorithm  was  ji  =  500  for  (pU2  and  105  for  </>Ul 
and  (f)us  • 

Figure  3.10  shows  the  rms  velocities  when  the  rms  component  is  included  in  the  objective 
function.  As  illustrated  in  Figure  3.10,  the  prediction  of  the  rms  velocities  improves  when  the 
transpiration  velocity  control  is  added;  however,  the  streamwise  rms  velocity  is  still  over-predicted 
near  the  wall.  Not  shown  for  this  simulation  is  the  mean  velocity  profile,  which  in  this  case  is  not 
as  good  as  the  mean  velocity  profile  that  is  achieved  in  the  previous  section  when  only  the  mean 
velocity  profile  is  targeted  by  the  controls.  If  shown,  it  would  lie  between  the  two  mean  velocity 
profiles  in  Figure  3.2  corresponding  to  control  by  wall  stress  only  and  control  by  wall  stress  plus 
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Figure  3.10:  Root-mean-square  of  velocity  fluctuations  with  objective  to  control  mean  flow  and 

rms  velocities.  - ,  control  includes  wall  stresses  and  transpiration; - ,  control  includes  wall 

stresses  only; .  ,  reference  profiles  from  [24]. 


transpiration  velocity.  Furthermore,  the  region  in  which  the  improved  predictions  occur  is  limited 
to  approximately  the  first  three  grid  cells  adjacent  to  the  wall. 

The  results  of  this  simulation  show  that  the  prediction  of  velocity  fluctuation  intensities  can  be 
improved  by  the  addition  of  a  wall-normal  velocity  approximate  boundary  condition.  But,  the  fact 
that  the  mean  velocity  profile  is  not  as  well  predicted  when  the  velocity  fluctuations  are  targeted 
through  the  objective  function  suggests  that  the  objectives  of  getting  the  correct  mean  velocity 
profile  and  the  correct  rms  velocities  may  be  competing  objectives. 

3.5  A  simpler  wall  model  from  LSE  of  wall  stress  data 

While  the  sub-optimal  control  strategy  for  generating  wall  stresses  could  be  used  as  a  wall  model 
for  coarse-grid  LES,  its  cost  is  approximately  20  times  greater  than  of  an  LES  on  the  same  grid 
compared  to  an  explicit  wall  stress  model  such  as  Eq.  (1.1).  Furthermore,  a  target  mean  velocity 
profile  must  be  provided  to  define  the  objective  function.  It  may  be  possible  to  lower  the  cost 
of  control  strategy,  but  that  possibility  is  not  investigated  here.  The  real  strength  of  the  optimal 
control  strategy  is  that  it  yields  wall  stress  boundary  conditions  that  work  for  coarse-grid,  high 
Reynolds  number  LES.  Thus,  a  reference  data  set  can  be  generated  against  which  new  wall  models 
can  be  compared.  The  reference  data  can  even  be  used  to  derive  new  wall  models.  Such  an  approach 
is  described  in  this  section. 

The  most  desirable  wall  stress  model  would  be  similar  to  Eq.  (1.1)  in  that  it  would  be  an  explicit 
relation  between  the  wall  stresses  and  the  velocity  field.  One  could,  for  instance,  require  the  wall 
stress  model  to  be  the  best  possible  mean  square  estimate  of  the  sub-optimal  wall  stress  as  a 
function  of  the  velocity  field  in  a  neighborhood  of  the  point  where  the  wall  stress  is  required.  This 
is  the  conditional  average  of  the  wall  stress  given  the  local  velocity  field  (a  conditional  average  is 
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necessary  because  the  wall  stress  may  have  a  stochastic,  or  unpredictable,  component  with  respect 
to  the  local  velocities).  It  is  denoted  by  z)|E),  where  E  is  a  vector  of  events.  In  the  present 

study,  E  will  be  a  vector  containing  the  local  velocity  field,  but  it  could  easily  contain  pressure, 
velocity  gradients,  quadratic  products,  or  any  other  quantities  which  might  characterize  the  wall 
stresses.  The  conditional  average  embodies  so  much  statistical  information  that  it  is  unlikely  that 
it  could  be  found  exactly,  but  it  can  be  approximated  by  its  linear  stochastic  estimate  (LSE),  given 
by  (see  [3]  for  instance): 

(Tg{x,z)\E)*f3{x,z)=LijEj  *  =  1,3,  i  =  l,2,3,...,JV,  (3.28) 

where  N  is  the  number  of  events  being  considered,  and  Lij  are  estimation  coefficients  relating 
to  Ej .  By  the  statistical  orthogonality  principle  [34],  the  mean  square  error  between  and  f$  is 
minimized  when  the  event  data  are  uncorrelated  with  the  error  e*  —  —  f$: 

(eiEk)  =  ((r%-f%)Ek)  =  0.  (3.29) 

Substituting  Eq.  (3.28),  the  estimation  coefficients  are  governed  by: 

(r%Ek)=Lij(EjEk).  (3.30) 

The  matrix  ( EjEk )  is  invertible  provided  the  events  are  linearly  independent.  Thus,  to  obtain  the 
LSE,  the  correlations  {r^Ek)  and  ( EjEk )  must  be  found  and  the  events  must  be  selected  that  best 
characterize  the  wall  stress. 

Though  the  technique  employed  here  is  essentially  the  same  as  that  of  Bagwell  et  al.  [7],  the 
results  are  different  since  the  reference  data  used  here  is  already  known  to  work  well  for  a  coarse  grid 
LES,  whereas  Bagwell’s  reference  data  comes  from  a  direct  numerical  simulation  at  low  Reynolds 
number.  As  will  be  shown  below,  an  event  field  consisting  of  the  nearby  velocities  is  sufficient  to 
yield  wall  models  of  the  form  (3.28)  that  have  greater  than  80%  correlation  with  the  optimal  wall 
stresses.  Moreover,  the  new  wall  models,  when  used  in  an  LES,  will  be  shown  to  reproduce  the 
results  of  the  sub-optimal  control  LES. 

3.5.1  LSE  predictions 

To  implement  the  LSE,  eight  velocity  fields,  well  separated  in  time,  and  their  sub-optimal  wall 
stresses  from  the  Rer  =  4000  simulation  discussed  in  Section  3.3  were  used.  In  implementing  the 
LSE,  it  was  found  that  better  results  were  obtained  by  estimating  only  the  fluctuating  part  of  the 
wall  stresses  from  the  fluctuating  part  of  the  velocity  field.  To  this  end,  the  instantaneous  plane 
averaged  wall  stress  was  subtracted  from  the  wall  stresses  and  the  average  over  the  first  off-wall 
plane  of  the  horizontal  velocity  was  subtracted  from  its  velocity  component  for  each  sample.  The 
LSE  can  be  written  as  a  convolution  sum  in  the  wall-parallel,  homogeneous  directions.  The  exact 
form  of  the  LSE  wall  model  for  the  fluctuating  part  of  the  wall  stresses  is  given  by: 

Tlz  fly  Tlx 

fu\m,n=  £  £  £  L\jk[Um-ij,n-k  “(«(•, |ft,*)>] 

k——nz  j— 1  Z— — Tlx 

7lZ+l  fly  flx  —  1 

X)  ]X  IX  ^ijkVJrn-ijtn-ky  (3.31) 

k=—nz  j— 1  X— — Tix 
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nz—  1  Tiy  nx+l 

^32  I m,n  =  X^  X^  Xy  ^ijk[um-i,j,n-k  ~  {u(',yi,  •))] 

k——nz  j= 1  i——nx 
riz  ny  nx 

+  Y1  X]  Lijkwm-i,j,n-k, 

k=~nz  j= 1  i=—nx 


(3.32) 


where  the  parameters  nx,ny,nz  determine  the  number  of  velocity  points  used  in  the  convolution 
sums  to  estimate  the  wall  stress  at  each  wall  location  (denoted  by  the  subscript  pair  m,  n).  Note  that 
the  summation  relating  v2  and  the  spanwise  velocity  w  has  different  indices  than  the  summation 
for  u  in  Eq.  (3.31).  This  is  due  to  the  staggered  grid.  Also  note  that  the  wall-normal  velocity  does 
not  appear  in  the  LSE  wall  model.  This  is  because  the  wall-normal  velocity  is  linearly  dependent 
on  the  wall-parallel  velocities  and  thus  cannot  be  used  as  an  independent  event  in  the  LSE.  The 
coefficients  Ln,L13,L31,  and  L33,  are  determined  by  requiring  that  the  error  be  orthogonal  to  the 
events  (velocities)  as  in  Eq.  (3.29). 

For  an  a  priori  comparison  of  the  wall  stress  fluctuations  predicted  by  the  LSE  and  those  of  the 
optimal  strategy  one  possibility  is  to  compute  the  correlation  coefficient: 


(Min 


w'\ 
i2  ) 


((fg)W(rg)2y^ 


(3.33) 


where  i  =  1, 3  and  the  average  denoted  by,  (•),  is  taken  over  all  of  the  samples.  Another  frequently 
used  quantity  for  comparison  is  the  relative  mean  square  error: 


Ri2  = 


((rg)2) 


(3.34) 


where,  again,  i  =  1,3. 

Table  3.1  shows  the  correlation  coefficients  and  relative  mean  square  error  for  several  choices  of 
the  parameters  nx,ny,  and  nz.  The  first  noteworthy  observation  is  that  the  correlations  between 
the  sub-optimal  wall  stress  fluctuations  and  those  predicted  by  the  LSE,  Eqs.  (3.31)  and  (3.32) 
increase  significantly  when  velocity  data  is  included  from  the  first  two  wall-parallel  planes  (compare 
Tiy  =  1  and  ny  =  2).  The  estimates  with  nx  =  15  and  nz  =  15  use  nearly  all  of  the  velocity  data 
in  each  wall-parallel  plane,  but  note  that  this  results  in  little  improvement  over  the  comparable 
nx  =  4,nz  =  4  results.  Figures  3.11-3.14  show  contour  plots  of  the  LSE  coefficients  for  the  case 
nx  =  2,  ny  =  2,  nz  =  2.  Each  figure  contains  two  plots,  the  top  plot  corresponding  to  the  first  plane 
of  velocity  data  and  the  bottom  plot  to  the  second.  The  plots  indicate  several  symmetries  in  the 
LSE  coefficients.  For  example  L11  is  symmetric  with  respect  to  reflection  in  the  spanwise  direction 
whereas  L13  is  anti-symmetric.  These  symmetries  could  be  exploited  to  reduce  the  number  of  free 
coefficients  in  a  wall  model,  but  we  have  not  explored  that  here. 

Bagwell  et  al  [7]  used  entire  planes  of  velocity  data  to  form  the  LSE,  but  as  results  in  Table  3.1 
suggest,  a  small  number  of  data  points  may  be  sufficient.  These  a  priori  results  suggest  that  the 
LSE  can  reproduce  the  sub-optimal  wall  stress  fluctuations  reasonably  well,  but  even  correlations 
in  excess  of  80%  do  not  guarantee  that  the  LSE  can  perform  well  as  a  wall  model.  This  can  only 
be  checked  via  an  actual  simulation. 


3.5.2  LSE  wall  model  results 

The  LSE  was  used  to  model  the  fluctuating  portion  of  the  wall  stresses  in  terms  of  the  velocity  field 
and  must  be  combined  with  a  model  for  the  mean  wall  stress  to  be  used  as  an  LES  wall  model. 
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Figure  3.11:  Contours  of  LSE  coefficient  Ln  for  nx  =  2 ,ny  =  2,  and  nz  =  2,  in  the  zz-plane.  The 
top  plot  is  for  the  first  plane  (j  =  1)  and  the  bottom  plot  is  for  the  second  plane. 
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Figure  3.12:  Contours  of  LSE  coefficient  L13.  Refer  to  Figure  3.11  for  details. 
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Figure  3.13:  Contours  of  LSE  coefficient  L31.  Refer  to  Figure  3.11  for  details. 
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Figure  3.14:  Contours  of  LSE  coefficients33.  Refer  to  Figure  3.11  for  details. 
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Table  3.1:  Correlation  coefficients  and  relative  mean  square  errors  for  the  LSE  wall  model  fluctu¬ 
ations  versus  the  sub-optimal  wall  stress  fluctuations  for  different  numbers  of  velocity  data  use  in 
the  estimates.  The  fifth  line,  set  apart  by  horizontal  lines,  is  for  the  coefficients  used  as  an  actual 
wall  model  in  Section  3.5.2. 
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The  actual  model  used  is: 


rf2\m,n=  <rr2)+f#|  m,„ 

732 1 ra, n  —  T32  |m,ra? 


(3.35) 

(3.36) 


where  ff2  and  are  given  by  Eqs.  (3.31)  and  (3.32),  respectively.  The  mean  wall  stress  (t^)  is 
found  by  assuming  that  the  plane-averaged  streamwise  velocity  at  the  first  off- wall  grid  location 
and  the  mean  wall  stress  are  related  by  the  logarithmic  law  of  the  wall: 


=  (Ti2>1/2 


2.41  log( 


)  +  5.2 


(3.37) 


By  using  the  LSE  wall  model  (3.35)  in  the  same  flow  discussed  in  Section  3.2.1,  with  Rer  =  4000 
and  the  same  resolution,  it  was  found  that  the  wall  models  based  on  the  one  plane  LSE,  ny  =  1 
resulted  in  mean  velocity  profiles  (not  shown  here)  that  were  not  as  good  as  the  profiles  from  the 
sub-optimal  simulation.  However,  the  wall  model  based  on  the  two  plane  LSE  was  able  to  reproduce 
nearly  exactly  the  results  of  the  sub-optimal  simulation.  Shown  in  Figure  3.15  is  the  mean  velocity 
profile  from  a  simulation  with  an  LSE  wall  model  of  the  form  (3.35)  with  nx  =  2,  ny  =  2,  and  nz  —  2. 
Simulations  with  LSE  models  based  on  smaller  stencils  of  velocity  data  did  not  work  as  well,  whereas 
LSE  models  based  on  larger  stencils  reproduced  the  sub-optimal  results.  The  turbulent  fluctuations 
are  also  in  perfect  agreement  with  those  of  the  sub-optimal  wall  stress  simulation;  see  Figure  3.16. 

These  results  are  remarkable  in  that  the  wall  stress  model  (3.35)  reproduces  the  results  of  the 
sub-optimal  simulation  at  a  cost  only  slightly  higher  than  that  of  a  simulation  with  no  wall  model. 
However,  the  model  given  by  Eq.  (3.35)  will  not  be  of  much  use  if  new  coefficients  need  to  be  derived 
for  different  Reynolds  numbers  or  for  different  grid  resolutions.  Fortunately,  we  have  found  that 
simulations  with  Eq.  (3.35)  perform  well  over  a  large  range  of  Reynolds  numbers  using  the  same 
coefficients  derived  from  the  Rer  =  4000  sub-optimal  data.  Figure  3.17  shows  the  results  of  using 
Eq.  (3.35)  at  Reynolds  numbers:  Rer  =  640, 4000,  and  20000.  The  simulations  were  conducted  on 
the  same  uniform  grid  as  the  ReT  =  4000  simulation  with  the  log  law,  Eq.  (3.37),  used  to  determine 
the  mean  wall  stress.  In  the  ReT  =  640  LES,  the  log  law  was  used  to  relate  the  plane-averaged 
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Figure  3.15:  Mean  velocity  profiles  for  Rex  =  4000.  - ,  sub-optimal  computation; 

Shifted  model  of  Piomelli  et  al.  [36]; - ,  LSE  model; - ,  {u)+  =  2.41  In  y+  +  5.2. 


Figure  3.16:  Root-mean-square  of  velocity  fluctuations  for  ReT  =  4000  and  uniform  32  x  33  x  32 
grid. - ,  sub-optimal  computation; - ,  LSE  model. 
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Figure  3.17:  Mean  velocity  profiles  from  LES’s  with  the  LSE  model  on  uniform  32  x  33  x  32  grid: 
- ,  ReT  =  640; - ,  Rer  =  4000; - ,  Rer  =  20000; - ,  (u)+  =  2.41  In  y+  +  5.2. 


streamwise  velocity  from  the  second  plane  to  the  mean  wall  stress  because  the  first  plane  is  in  the 
buffer  region  at  y+  «  20. 

Similarly,  the  same  coefficients  were  used  with  Eq.  (3.35)  in  a  simulation  with  a  refined  grid  at 
Rer  =  20000.  The  number  of  cells  was  doubled  in  each  direction  to  64  x  65  x  64.  The  resulting 
mean  flow  profile  is  shown  in  Figure  3.18.  The  log  region  intercept  is  still  slightly  underpredicted, 
but  the  mean  flow  now  exhibits  a  wake-like  behavior  near  the  channel  center  as  has  been  observed 
in  high  Reynolds  number  channel  simulations  in  which  the  near- wall  region  is  resolved  [24], 

Unfortunately,  this  simple  linear  model  is  not  going  to  be  a  panacea.  Figure  3.19  shows  the 
mean  velocity  profiles  for  several  channel  flow  LES’s  at  ReT  =  4000,  all  using  the  same  number 
of  grid  points  as  the  simulations  discussed  above  and  using  the  simple  linear  wall  stress  model 
previously  derived.  In  each  case  some  reasonable  modification  has  been  made.  For  instance,  a  fully 
conservative  fourth  order  finite  difference  scheme  was  used,  and,  as  shown  in  the  figure,  the  mean- 
velocity  is  under-predicted.  To  test  the  effects  of  the  numerics  on  the  efficacy  of  the  wall  model, 
two  different  things  were  tried:  stretching  the  grid  in  the  wall-normal  direction  and  modifying  the 
dynamic  procedure  as  suggested  by  [13].  As  Figure  3.19  shows,  the  simple  linear  wall  stress  model 
performs  worse  in  every  one  of  these  cases  than  in  the  original  simulation  for  which  it  was  designed. 

3.6  Conclusion 

A  suboptimal  control  strategy  has  been  successfully  applied  to  a  coarse  grid  LES  of  a  channel 
flow  at  high  Reynolds  number.  The  two  objectives  of  this  study  are:  1)  to  demonstrate  that  a 
control  strategy  can  determine  approximate  wall  boundary  conditions  that  result  in  an  accurate 
LES,  and  2)  to  find  a  simple  wall  model  using  these  results.  In  both  cases,  the  work  can  be  judged 
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Figure  3.18:  Mean  velocity  profiles  from  LES’s  for  Rer  =  20000  and  different  uniform  grids. - , 

64  x  65  x  64  grid  points  and  LSE  model; - —  ,  32  x  33  x  32  grid  points  and  LSE  model;  — i —  , 

64  x  65  x  64  grid  points  and  shifted  model  of  Piomelli  et  al.  [36]; - ,  (u)+  =  2.41  lny+  +  5.2. 


Figure  3.19:  Mean  velocity  profiles  using  fixed,  simple  linear  model  for  the  wall  stresses.  .  , 

logarithmic  reference  profile  =  2.41  In  y+  +  5.2;  - ,  model  reproduces  mean  profile  when 

used  in  same  setting  that  it  was  derived; - ,  same  model  with  fourth-order  finite  differences; 

- ,  same  model  with  modified  dynamic  procedure  as  in  [13];  -  o  -  ,  same  model  with  stretched 

wall-normal  grid. 
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a  partial  success.  The  sub-optimally  controlled  LES,  using  both  wall  stresses  only  and  wall  stresses 
and  transpiration  velocity,  is  able  to  generate  mean  velocity  profiles  in  good  agreement  with  the 
logarithmic  law.  It  is  also  determined  that  the  wall  stresses  provide  most  of  the  control  authority  in 
this  situation.  However,  the  rms  velocity  fluctuations,  particularly  the  streamwise  component  near 
the  wall,  are  not  well  predicted  using  this  scheme  when  the  objective  function  includes  the  mean 
flow  only.  When  the  objective  function  is  enhanced  to  include  a  component  targeting  the  reference 
rms  velocity  fluctuations  of  Kravchenko  [24],  modest  improvement  is  observed,  but  the  streamwise 
rms  velocity  remains  overpredicted  near  the  wall.  Furthermore,  the  objective  of  matching  the  rms 
velocity  fluctuations  appears  to  compete  with  the  objective  to  match  the  mean  velocity,  resulting 
in  a  slightly  less  accurate  mean  profile. 

The  data  generated  from  the  simulations  using  wall  stress  control  only  have  been  used  to  derive 
a  LSE  model  which  predicts  the  wall  stresses  from  the  near-wall  velocity  field.  The  LSE  model 
is  found  to  work  well  over  a  wide  range  of  Reynolds  numbers  using  the  same  grid  and  numerical 
methods.  However,  when  applied  to  LES  at  the  same  Reynolds  number  that  the  training  data  are 
generated,  but  with  different  numerical  methods,  the  same  model  is  found  to  perform  poorly.  This 
provides  additional  evidence  that  non-physical  errors,  namely  numerical  and  SGS  modeling  errors, 
play  a  critical  role  in  the  wall  modeling  problem  and  must  be  addressed  for  an  adequate  solution 
to  be  found. 
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Chapter  4 

New  Directions  in  Control-Based  Wall 
Modeling 


Many  important  lessons  were  learned  from  the  control  based  wall  modeling  work  discussed  in  the 
previous  chapter.  Unfortunately,  this  approach  is  impractical  due  to  the  high  computational  cost 
required  for  the  suboptimal  control,  which  requires  both  the  solution  of  adjoint  equations  and 
many  iterations  to  achieve  convergence  in  the  wall  stresses.  Furthermore,  the  cost  function  is  based 
on  known  target  data,  making  the  model  non-predictive.  The  LSE  models  generated  from  such 
computations  are  too  sensitive  to  the  numerical  parameters  to  construct  a  universal  LSE  coefficient 
database.  Thus,  a  low-cost,  robust  wall  model  is  needed  to  achieve  the  accuracy  of  the  sub-optimal 
control  technique  without  an  a  priori  target  solution.  To  this  end,  a  cost  function  based  on  a 
RANS  solution  will  be  constructed  in  Section  4.1  to  make  the  model  predictive,  and  in  Section  4.2, 
the  problem  will  be  formulated  in  an  optimal  shape  design  setting  in  an  attempt  to  reduce  the 
computational  cost.  Some  test  results  and  discussions  are  presented  in  Sections  4.3  and  4.4. 


4.1  Cost  function 


In  order  to  make  the  wall  model  predictive,  an  easy-to-evaluate  cost  function  near  the  wall  using 
quantities  not  known  a  priori  must  be  defined.  To  this  end,  a  RANS  model  is  used  to  provide  the 
target  velocity.  This  is  motivated  by  the  recognition  that  the  near- wall  region  of  a  high  Reynolds 
number  boundary  layer  is  more  appropriately  modeled  by  RANS  than  by  a  coarse  grid  LES  with 
filter  length  larger  than  the  integral  scale  of  the  turbulence. 

In  the  present  work,  the  RANS  model  is  obtained  from  the  simplified  version  of  the  TBL  equation 
model  introduced  in  Chapter  2, 


A 

dy  L 


+  vt  (»)) 


dui 

dy  J 


1  dp 

p  dx<i 
2 


1,3 


LES 


t'i  ( y )  =  Kvy+  (l  -  e  y+!A^j  ,  y+  =  yuT/ u. 


(4.1) 


These  equations  model  all  Reynolds  stresses  through  a  damped  mixing-length  eddy  viscosity,  and 
explicitly  account  for  the  pressure  gradient  which  is  assumed  constant  across  the  wall  layer  and  is 
imposed  by  the  LES.  To  complete  the  model,  a  no  slip  condition  is  applied  at  the  wall  and  the 
outer  boundary  is  set  to  be  the  LES  velocity.  The  resulting  velocity  profile  should  be  interpreted 
as  the  ensemble  averaged  velocity  profile  given  the  local  LES  state.  It  can  therefore  be  expected 
that,  on  average,  the  resolved  LES  should  match  the  RANS  solution  near  the  wall.  Note  that  this 


35 


model  is  chosen  for  simplicity  in  this  initial  attempt,  and  there  are  likely  better  models  for  this 
application  that  will  be  explored  in  future  work. 

In  an  overlapped  region  consisting  of  N  LES  grid  points  in  the  wall-normal  direction,  cost 
functions  are  devised  to  match  the  LES  and  RANS  solutions  on  average.  An  attractive  method  in 
a  statistically  stationary  flow  would  be  to  use  a  running  time  average  to  provide  the  target  velocities. 
However,  if  the  control  authority  is  restricted  to  the  current  time,  this  approach  becomes  impractical 
since  the  flow  at  the  current  time  would  contribute  only  a  small  fraction  of  the  total  cost  function. 
This  makes  it  difficult  to  determine  the  control  since  the  cost  function  is  insensitive  to  it.  If  the 
control  is  explicitly  computed  as  a  function  of  time,  then  adjoint  equations  have  to  be  integrated 
backward  in  time  to  find  a  correct  solution  over  a  sufficiently  large  time  window  which  contains 
enough  statistical  samples. 

An  alternative  is  to  use  the  current  state  as  the  statistical  sample.  Thus,  the  first  cost  function 
is  defined  to  be  the  £2  difference  between  the  LES  and  RANS  states: 


r  N 

Jc2  =  JSY1  ((^RANS.llyn  -  «LES,l|j/„)2  +  (URANS,3|»„  -  ULES^IyJ2)  dS, 


(4.2) 


where  S  is  the  surface  and  yn  are  the  locations  of  the  n  overlap  points.  In  this  way,  a  sufficient 
number  of  samples  of  the  flow  state  are  used  to  make  a  meaningful  average.  Also,  the  cost  function 
is  based  only  on  quantities  at  the  current  time  step,  so  no  history  information  is  required.  This 
type  of  cost  function  is  also  compatible  with  the  gradient  evaluation  methods  used  in  this  work 
(see  Section  4.2). 

Other  cost  functions  can  also  be  formulated  for  this  problem.  A  cost  function  based  on  the 
average  deviation  of  the  LES  and  RANS  is: 


JA  =  ^  ((uRANS,l|yn  -  ULES,l|</J  +  (^RANS,3|yn 


j 


This  cost  function  is  similar  to  that  used  in  [33].  However,  as  shown  in  Section  4.3,  this  cost 
function  performs  quite  poorly.  Analysis  of  its  gradients  indicates  that  they  do  not  capture  the 
sign  information  correctly  in  some  regions  (gradient  computation  will  be  discussed  in  the  next 
section).  In  order  to  retain  more  information  and  move  in  the  direction  of  feedback  control,  a 
signed  cost  function  has  also  been  used: 


Js 


((uRANS,l|y„ 


-  «LES,l|y»)  +  (uRANS,3|y„ 


«LES,3|yJ)  dS. 


(4.4) 


When  this  cost  function  is  used,  the  control  strategy  is  shifted  to  force  the  cost  function  to  zero 
rather  than  minimizing  it.  It  was  thought  that  this  approach  might  better  take  advantage  of  the 
method  being  used  for  gradient  evaluation,  but  it  only  resulted  in  a  moderate  improvement  (see 
Section  4.3). 

The  choice  of  N  in  Eqs.  (4.2)-(4.4)  should  be  made  to  include  as  many  matching  layers  as 
possible  while  remaining  in  the  region  where  the  RANS  model  is  a  reasonable  approximation  for 
the  given  local  flow.  Furthermore,  the  LES  velocity  too  close  to  the  wall  may  involve  large  errors 
(see  [12])  and  thus  is  not  suitable  as  a  RANS  boundary  condition.  In  the  calculations  presented  in 
this  article,  N  has  been  chosen  to  be  three. 

Two  important  points  should  now  be  noted.  First,  while  all  the  cost  functions  here  are  based 
on  matching  RANS  and  LES  velocities,  other  quantities  could  also  be  used.  These  could  include 
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matching  vorticity  or  energy  fluxes  with  suitable  models.  Second,  it  may  not  be  possible  or  desirable 
to  reduce  the  cost  function  to  zero.  Doing  so  could  artificially  reduce  the  turbulence  fluctuations 
of  the  flow.  Also,  if  an  inexpensive  scheme  is  required,  it  may  not  be  possible  to  fully  optimize  the 
solution.  Thus,  the  cost  function  must  act  as  a  suitable  quantity  for  feedback  regulation,  rather 
than  for  minimization. 


4.2  Optimization  using  shape  design  techniques 

Optimal  shape  design  consists  of  a  set  of  techniques  for  optimizing  a  shape  to  achieve  an  engineering 
objective  (e.g.  [31]).  Several  approaches  have  been  developed  in  this  field  that  have  had  some  success 
in  reducing  the  computational  expense  of  the  optimization  procedure.  In  an  attempt  to  bring  these 
techniques  to  bear,  the  wall  modeling  problem  is  formulated  in  this  framework. 

In  general,  the  formulation  is  to  consider  a  partial  differential  equation  A  (U,  q,a)  —  0  in  a  region 
fi  satisfying  boundary  conditions  b  ({7,  q,  a)  =0  on  90.  The  optimization  is  performed  to  determine 


min  { J  (£/,  q,  a)  :  A  (£7,  q,  a)  =  0  Vx  E  0,  b  ( U ,  q,  a)  ~  0  Vx  €  90}  (4.5) 

for  some  cost  function  J  (£7,  g,  a).  In  this  formulation,  U  is  the  state,  q  the  shape,  and  a  are  the 
control  variables.  The  gradient  of  the  cost  function  with  respect  to  the  control  variables  is  then: 

da  da  dq  da  ^  dU  dq  da 

The  standard  technique  for  solving  this  equation  is  to  use  an  adjoint  method  interfaced  with  a 
gradient  minimization  technique.  But,  as  previously  noted,  this  can  be  expensive  and  present  data 
storage  difficulties  in  time-accurate  computations.  Since  it  is  the  last  term  in  Eq.  (4.6)  that  requires 
the  adjoint  evaluation,  Mohammadi  &  Pironneau  [31]  suggest  the  following  assumption  when  the 
controls  and  the  cost  function  share  the  same  support: 


dJ  dJ  dJ  dq 

_  _ |_ _ ± 

da  da  dq  da 


(4.7) 


This  assumption  is  called  the  method  of  incomplete  sensitivities  since  the  sensitivity  to  the  state 
gradient  is  ignored.  The  use  of  this  method  has  been  explored  in  this  work  since  it  has  produced 
positive  results  in  the  optimization  of  aerodynamic  shapes.  For  examples,  see  [29,  30,  31],  although 
these  are  all  steady,  two-dimensional  applications.  Since  no  rigorous  proof  on  the  applicability  of 
this  technique  exists  and  its  usefulness  is  based  on  purely  empirical  studies,  it  was  not  known  how 
well  it  would  perform  in  a  full  LES.  Furthermore,  the  present  cost  function  is  not  defined  exactly 
on  the  support  of  the  control,  although  it  is  defined  in  a  small  neighborhood  of  the  control.  While 
these  factors  will  produce  errors,  the  gradient  evaluation  needs  only  accurately  predict  the  sign  of 
the  gradient  and  capture  to  some  degree  the  difference  in  magnitudes  of  the  derivatives  with  respect 
to  different  control  parameters.  A  goal  of  this  work  is  to  determine  if  the  amount  of  information 
contained  in  this  gradient  is  sufficient  for  application  to  wall  boundary  conditions. 

In  order  to  apply  the  incomplete  sensitivity  assumption,  the  control  must  be  related  to  shape 
design  parameters.  B-splines  spaced  evenly  along  the  surface  (although  not  enough  to  form  a  com¬ 
plete  basis)  are  used  to  parameterize  deformations  normal  to  the  surface.  The  control  parameters, 
a*,  are  then  the  spline  amplitudes.  The  gradient  of  the  cost  function  with  respect  to  these  param¬ 
eters  can  be  computed  using  finite  differences  by  perturbing  each  parameter  by  a  small  value,  e, 
and  then  using  Eq.  (4.7)  to  evaluate  the  gradient  based  on  the  current  state  information.  It  is  not 
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necessary  to  recompute  the  actual  geometry  or  grid  because  all  the  state  variables  of  interest  can 
be  stored  and  matched  to  the  new  surface.  The  parameter  e  is  chosen  a  priori  by  making  it  small 
enough  such  that  the  gradient  values  are  independent  of  it. 

Once  the  cost  function  gradient  is  known,  the  new  spline  amplitudes  can  be  computed  by 


„*+i 


al 


dJ 
Pdai ’ 


(4.8) 


where  p  is  a  descent  parameter  set  in  advance  and  k  is  the  iteration  count.  The  new  shape  is 
computed  by  adding  the  surface  perturbations  to  the  previous  shape.  To  relate  this  to  the  wall 
stresses,  the  RANS  model  is  used  to  compute  the  correction  to  the  equivalent  slip  velocity  on  the 
original  surface: 

=  /RANS,i  {l/new )  5  i  =  1,  3,  (4*9) 

where  /  stands  for  the  RANS  model  given  by  Eq.  (4.1).  This  approach  is  inspired  by  a  Taylor 
series  expansion  about  the  wall  [31].  In  this  way,  it  is  not  necessary  to  change  the  computational 
geometry  of  the  LES. 

The  total  slip  velocity  is  given  by  adding  the  correction  ucw  i  to  the  old  wall  slip  velocity.  Cor¬ 
rected  wall  stresses  can  then  be  computed  directly  by  definition 


,i  =  + 


-tc 


Re  &X2  ’ 


(4.10) 


where  Ax2  is  the  local  wall  normal  grid  spacing. 

While  this  approach  avoids  the  evaluation  of  a  set  of  adjoint  equations,  iterations  are  still 
required  to  converge  the  solution.  Additional  function  evaluations  are  also  often  used  to  determine 
an  optimal  choice  for  p  at  each  iteration.  In  order  to  make  the  wall  model  practical,  these  costs 
must  be  avoided.  Therefore,  no  iterations  are  performed  at  each  time  step.  The  cost  function 
gradients  are  computed  and  used  in  a  feedback  manner  to  provide  a  correction.  Every  a*  is  reset 
to  zero  at  each  time  step.  Also,  p  is  taken  to  be  a  fixed  parameter  similar  to  the  gain  in  a  feedback 
controller.  To  make  up  for  some  of  this  lost  information,  a  predictor-corrector  approach  to  the 
control  algorithm  will  be  used.  This  is  done  by  using  Eq.  (4.1)  to  compute  a  prediction  of  the 
wall  stress  before  the  optimization  is  used.  It  is  expected  that  the  prediction  will  account  for  the 
missing  physics  in  the  coarse  grid  LES  while  the  optimization  will  correct  for  the  numerical  and 
SGS  modeling  errors.  While  this  approach  must  be  classified  as  sub-optimal,  it  is  still  reasonable 
to  expect  a  cost  function  reduction  if  at  each  time  step  the  LES  velocity  is  forced  in  the  direction 
of  the  reduced  cost  function. 


4.3  Results 

The  application  of  this  method  to  the  trailing  edge  flow  simulated  previously  by  Wang  &  Moin 
[41,  42]  has  produced  mixed  results.  The  first  goal  is  to  justify  the  incomplete  sensitivities  assump¬ 
tion.  The  £2  cost  function  history  is  shown  in  Figure  4.1.  While  the  average  value  is  reduced 
approximately  15%  from  the  initial  value,  this  is  not  completely  out  of  the  range  of  the  cost  func¬ 
tion  fluctuations.  It  is  therefore  inconclusive  regarding  the  validity  of  the  assumption.  As  shown 
in  Figure  4.2,  the  predicted  wall  stress  matches  the  full  LES  wall  stress  quite  well  in  some  regions 
for  the  £2  and  signed  cost  functions,  but  performs  poorly  in  other  regions.  The  separation  point 
is  predicted  reasonably  accurately  for  both  these  cost  functions.  As  previously  indicated,  the  av¬ 
erage  cost  function  performed  more  poorly.  Figure  4.3  contains  a  comparison  between  the  £2  cost 
function  results  and  the  predictor  alone.  The  new  results  are  much  better  in  the  region  near  the 
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x/h 

Figure  4.3:  Time  averaged  skin  friction  over  the  airfoil  surface: - ,  £2  cost  function; 

predictor  only; - ,  full  LES  [41]; - ,  TBL  equation  model  [42]. 


skin  friction  peak,  although  they  produce  a  less  smooth  skin  friction  profile,  and  rather  large  errors 
remain  in  part  of  the  adverse  pressure  gradient  region.  Overall,  the  model  demonstrates  some 
improvement  over  the  simple  wall  model  used  as  a  predictor,  but  is  less  accurate  than  the  full  TBL 
equation  model  used  in  [42]. 

Comparison  of  the  velocities  between  the  full  LES  and  wall  modeled  LES  (based  on  the  £2  cost 
function,  which  produced  the  best  results)  are  quite  good.  As  shown  in  Figure  4.4,  the  coarse  grid 
LES  is  able  to  match  the  resolved  LES  very  closely.  The  main  (moderate)  discrepancy  occurs  in 
the  turbulent  intensities  near  the  wall.  This  is  not  unreasonable  since  these  quantities  were  not 
included  in  the  cost  function  and  it  may  in  fact  not  be  possible  to  capture  these  regions  accurately 
because  the  LES  grid  does  not  resolve  the  intensity  peak.  When  compared  to  the  results  of  [42] 
using  only  the  predictor,  the  results  are  found  to  be  comparable  and  in  fact  are  worse  for  the  two 
cost  functions  not  shown.  Therefore,  it  is  difficult  to  draw  definitive  conclusions  about  the  effect  of 
the  gradient  based  optimization  procedure  on  the  velocity  field. 


4.4  Channel  flow  analysis 

In  order  to  evaluate  the  proposed  wall  model  in  a  more  controlled  environment,  the  algorithm  has 
been  implemented  in  the  plane  channel  LES  described  in  Chapter  3.  This  is  a  simpler  and  well 
known  case,  so  the  model  can  be  more  readily  analyzed.  It  was  immediately  noticed  that,  unlike 
the  trailing  edge  case,  the  cost  function  gradients  could  not  be  made  independent  of  the  small 
parameter  e  used  in  the  finite-difference  computation.  The  gradients  monotonically  decreased  with 
e  until  they  reached  a  value  of  zero.  This  result  indicated  that  the  incomplete  sensitivity  approach 
did  not  accurately  capture  the  gradients  in  the  channel,  since  non- zero  gradients  were  observed  in 
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Figure  4.4:  Mean  velocity  magnitude  profiles  (upper  figure)  and  streamwise  component  of  turbu¬ 
lence  intensities  (lower  figure)  at  several  trailing  edge  stations: - ,  £2  cost  function; - ,  full 

LES  [41],  Locations  are  those  indicated  in  Figure  4.2.  T.E.  is  the  trailing  edge  point. 
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the  sub-optimally  controlled  channel.  The  following  analysis  is  used  to  explain  these  results,  as 
well  as  the  difficulties  encountered  with  this  method  in  the  trailing  edge  geometry. 

Consider  a  cost  function  of  form 

J(a)  =  f  f(u(a))dS.  (4.11) 

Js 

Since  in  the  current  framework,  the  shape  and  shape  deformations  are  defined  in  two  dimensions, 
the  surface  can  be  parameterized  by  taking  the  y  coordinates  as  a  function  of  rr,  i.e.  y  =  g(x). 
Then  the  cost  function  becomes 

J(a)  =  f  f{u(x-,a))yj  1  +  gf2(x)dx.  (4.12) 

Jo 

Consider  a  perturbation  to  this  surface  parameterized  by  eh(x).  In  the  current  context,  h(x)  would 
correspond  to  the  spline  and  e  to  the  small  change  in  the  control  parameter.  The  new  cost  function 
is  computed  by  considering  its  sensitivity  to  geometry  only,  so 

J(a  +  e)  =  [  f(u(x\a))\/ 1  +  (g'(x)  +  eh'(x)y*dx.  (4.13) 

Jo 

By  using  a  Taylor  series  expansion,  one  obtains  to  0(e) : 

y/l  +  (0'(x)  +  e/i'(a:))2  «  y/i  TVH?)  +  e(l  +  9l2{x))~1/2g'{x)h'(x).  (4.14) 

When  the  gradient  is  computed  by  taking  ( J(a  +  e)  —  J(a))/e,  the  resulting  term  is 

^  «  J  f{u(x-,a))(l  +  gl2(x))~1^g'(x)h'{x)dx.  (4.15) 

This  expression  explains  the  observed  cost  function  gradients.  First,  it  has  been  demonstrated 
in  both  the  trailing  edge  and  channel  flows  that  in  regions  where  the  surface  is  flat,  the  gradients 
are  zero.  This  is  clear  since  in  these  regions,  g*(x)  =  0.  A  similar  observation  occurs  in  areas  where 
the  surface  is  a  straight  line.  This  is  because  gf(x )  is  constant  and,  in  this  case,  h(x)  is  symmetric, 
meaning  that  whenever  h'(x)  >  0,  there  is  a  corresponding  x\  such  that  hf(x i)  =  —h*(x).  Thus, 
unless  f(u(x\a))  has  a  very  large  change  between  x  and  xi,  since  gf(x)hf(x)  +  gf (xi)hf  (xi)  =  0  the 
gradient  will  be  very  small. 

Finally,  it  has  been  observed  that  in  regions  of  curvature  away  from  the  direction  of  perturbation 
and  for  a  positive  definite  /(u(x;a))  (such  as  the  £2  cost  function),  the  gradient  is  always  positive. 
This  can  be  seen  by  examining  the  product  gt(x)h,(x).  In  these  regions,  gf(x)  is  always  negative  and 
increases  monotonically  in  magnitude.  By  the  symmetry  of  h(x),  the  regions  where  h'(x)  is  positive 
correspond  to  gf(x)  having  a  smaller  magnitude,  and  the  regions  where  hr(x)  is  negative  correspond 
to  g((x)  having  a  greater  magnitude.  Thus,  the  positive  contribution  is  greater  in  magnitude  than 
the  negative  contribution,  and  hence  the  gradient  is  positive  since  f(u(x;a))  is  positive  and  varies 
less  than  the  curvature. 

The  sensitivity  computed  by  this  method  is  then  almost  exclusively  dependent  on  the  curva¬ 
ture  of  the  function  whose  information  is  contained  in  g*(x).  It  is  difficult  to  determine  how  this 
information  could  be  useful  in  changing  the  state  u  such  that  the  given  cost  function  is  minimized 
in  a  rigorous  and  well  defined  manner.  For  any  cost  function  defined  as  above,  the  incomplete 
sensitivity  method  will  act  in  a  way  directly  related  to  the  curvature  of  the  surface.  If  a  correlation 
exists  between  reducing  this  curvature  and  reducing  the  cost  function,  the  method  may  produce 
reasonable  results.  However,  there  is  no  reason  to  believe  that,  in  general,  reducing  surface  cur¬ 
vature  will  be  helpful  in  wall  modeling.  In  fact,  as  experience  in  the  channel  has  demonstrated,  a 
region  of  no  curvature  still  requires  control  to  obtain  an  accurate  solution.  Therefore,  it  is  likely 
that  an  alternative  method  must  be  found  for  the  general  application  of  a  wall  model. 
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4.5  Conclusions  and  future  work 

Wall  modeling  using  control  theory  is  a  promising  new  approach  for  developing  robust  wall  models 
which  account  for  not  only  the  unresolved  flow  physics  but  also  numerical  and  SGS  modeling 
errors.  In  this  chapter,  a  methodology  has  been  proposed  to  overcome  the  deficiencies  of  the  model 
by  Nicoud  et  al  [33],  described  in  Chapter  3,  and  make  the  control-based  wall  model  predictive  and 
practical  in  terms  of  computational  expense.  Two  critical  components,  namely  the  use  of  RANS 
velocity  profiles  as  the  near-wall  LES  target  in  the  cost  function  and  the  incomplete  sensitivity 
method  for  gradient  evaluation  have  been  examined  and  tested  in  a  turbulent  trailing  edge  flow. 

Based  on  the  test  results,  it  is  clear  that  the  assumption  of  incomplete  sensitivities  is  not 
appropriate  for  LES  wall  models  with  the  type  of  cost  function  considered  in  this  work.  This  is  at 
least  partly  due  to  the  cost  function  measuring  the  LES  state  in  the  flow  and  not  at  the  wall.  A 
cost  function  that  is  more  sensitive  to  the  geometry  could  be  better  suited,  but  it  is  unclear  how 
to  formulate  such  a  cost  function  for  a  wall  model.  Furthermore,  there  is  evidence  suggesting  that 
in  applications  similar  to  this,  the  gradient  calculated  with  incomplete  sensitivities  may  have  not 
only  incorrect  magnitude  but  also  incorrect  sign  [26].  Clearly,  a  more  accurate  means  is  needed  to 
compute  the  gradient. 

The  use  of  a  cost  function  matching  a  RANS  profile  near  the  wall  may  however  prove  useful  in 
LES  wall  modeling.  It  has  a  solid  physical  basis,  although  the  RANS  model  used  here  is  rather 
rudimentary.  More  robust  RANS  models,  such  the  k-u  model  are  being  considered.  In  addition  to 
choosing  an  appropriate  RANS  model,  the  choice  of  matching  quantities  is  also  an  important  factor 
in  the  performance  of  the  model.  Matching  LES  and  RANS  velocities  may  prove  not  to  be  the  best 
quantity  to  minimize  for  optimal  performance  of  the  model.  Cost  functions  based  on  vorticity  or 
energy  could  better  account  for  dynamics  that  are  more  important  to  the  large  scales  in  the  LES. 
An  investigation  of  these  cost  functions  and  implementation  of  a  RANS  model  is  underway  in  a 
channel  flow. 


43 


Acknowledgments 


This  work  was  supported  by  the  Air  Force  Office  of  Scientific  Research  through  contract  number 
F49620-00- 1-0111  (Dr.  Thomas  Beutner,  program  manager).  Computer  time  was  provided  by 
NAS  at  NASA  Ames  Research  Center  and  the  DoD’s  High  Performance  Computing  Modernization 
Program  though  ARL/MSRC.  The  authors  would  like  to  thank  Professors  Javier  Jimenez,  Peter 
Bradshaw,  and  Bijan  Mohammadi  for  many  helpful  discussions. 


44 


Bibliography 


[1]  F.  Abergel  and  R.  Temam.  On  some  control  problems  in  fluid  mechanics.  Theor.  Comput. 
Fluid  Dyn 1:303-325,  1990. 

[2]  E.  Achenbach.  Distribution  of  local  pressure  and  skin  friction  around  a  circular  cylinder  in 
cross-flow  up  to  Re  =  5  x  106.  J.  Fluid  Mech 34:625-639,  1968. 

[3]  R.J.  Adrian,  B.G.  Jones,  M.K.  Chung,  Y.  Hassan,  C.K.  Nithianandan,  and  A.  Tung.  Approx¬ 
imation  of  turbulent  conditional  averages  by  stochastic  estimation.  Phys.  Fluids  A,  1:992-998, 
1989. 

[4]  J.S.  Baggett,  J.  Jimenez,  and  A.G.  Kravchenko.  Resolution  requirements  in  large-eddy 
simulation  of  shear  flows.  Annual  Research  Briefs ,  Center  for  Turbulence  Research,  NASA 
Ames/Stanford  Univ.,  pp.  51-66,  1997. 

[5]  J.S.  Baggett,  F.  Nicoud,  B.  Mohammadi,  T.  Bewley,  J.  Gullbrand,  and  O.  Botella.  Sub-optimal 
control  based  wall  models  for  LES  -  including  transpiration  velocity.  Proc.  Summer  Program , 
Center  for  Turbulence  Research,  NASA  Ames/Stanford  Univ.,  pp.  331-342,  2000. 

[6]  T.G.  Bagwell.  Stochastic  Estimation  of  Near- Wall  Closure  in  Turbulence  Models.  PhD  thesis, 
University  of  Illinois  at  Urbana-Champaign,  1994. 

[7]  T.G.  Bagwell,  R.J.  Adrian,  R.D.  Moser,  and  J.  Kim.  Improved  approximation  of  wall  shear 
stress  boundary  conditions  for  large  eddy  simulation,  in  Near-Wall  Turbulent  Flows ,  edited  by 
R.M.C.  So,  C.G.  Speziale  and  B.E.  Launder,  Elsevier  Science,  New  York,  pp.  265-275,  1993. 

[8]  E.  Balaras,  C.  Benocci,  and  U.  Piomelli.  Two-layer  approximate  boundary  conditions  for 
large-eddy  simulations.  AIAA  Jounal ,  34:1111-1119,  1996. 

[9]  T.  Bewley,  H.  Choi,  R.  Temam,  and  P.  Moin.  Optimal  feedback  control  of  turbulent  channel 
flow.  Annual  Research  Briefs ,  Center  for  Turbulence  Research,  NASA  Ames/Stanford  Univ., 
pp.  3-14,  1993. 

[10]  T.R.  Bewley.  Flow  control:  new  challenges  for  a  new  renaissance.  Progress  in  Aerospace 
Sciences ,  37:21-58,  2001. 

[11]  M.  Breuer.  A  challenging  test  case  for  large  eddy  simulation:  high  Reynolds  number  circular 
cylinder  flow.  Int.  J.  Heat  Fluid  Flow ,  21:648-654,  2000. 

[12]  W.  Cabot.  Near-wall  models  in  large-eddy  simulations  of  flow  behind  a  backwards  facing 
step.  Annual  Research  Briefs ,  Center  for  Turbulence  Research,  NASA  Ames/Stanford  Univ., 
pp.  199-210,  1996. 


45 


[13]  W.  Cabot  and  P.  Mom.  Approximate  wall  boundary  conditions  in  the  large-eddy  simulation 
of  high  Reynolds  number  flow.  Flow  Turb.  Combust ,  63:269-291,  1999. 

[14]  Y.  Chang  and  S.S.  Collis.  Active  control  of  turbulent  channel  flows  based  on  large-eddy 
simulation.  ASME  Paper  No.  FEDSM-99-6929 ,  1999. 

[15]  H.  Choi,  R.  Temam,  P.  Moin,  and  J.  Kim.  Feedback  control  for  unsteady  flow  and  its  appli¬ 
cation  to  the  stochastic  burgers  equation.  J.  Fluid  Mech. ,  253:509-543,  1993. 

[16]  J.W.  Deardorff.  Numerical  study  of  three  dimensional  turbulent  channel  flow  at  large  Reynolds 
numbers.  J.  Fluid  Mech.,  41:453-480,  1970. 

[17]  M.  Germano.  A  dynamic  subgrid-scale  eddy  viscosity  model.  Phys.  Fluids  A,  3:1760-1765, 
1991. 

[18]  G.  Grotzbach.  Direct  numerical  and  large  eddy  simulation  of  turbulent  channel  flows,  in 
Encyclopedia  of  Fluid  Mechanics ,  edited  by  N.P.  Cheremisinoff,  Gulf  Publ.,  West  Orange,  NJ, 
pp.  1337-1391,  1987. 

[19]  F.H.  Harlow  and  J.E.  Welch.  Numerical  calculation  of  time-dependent  viscous  incompressible 
flow  of  fluid  with  free  surface.  Phys.  Fluids,  8:2182-2189,  1965. 

[20]  G.  Hoffman  and  C.  Benocci.  Approximate  wall  boundary  conditions  for  large  eddy  simulations, 
in  Advances  in  Turbulence  V,  edited  by  R.  Benzi,  Kluwer,  Berlin,  pp.  222-228,  1995. 

[21]  J.  Jimenez  and  R.D.  Moser.  LES:  where  we  are  and  what  we  can  expect.  AIAA  J.,  38:605-612, 

2000. 

[22]  J.  Kim,  P.  Moin,  and  R.  Moser.  Turbulence  statistics  in  fully  developed  channel  flow  at  low 
Reynolds  number.  J.  Fluid  Mech.,  177:133-166,  1987. 

[23]  A.G.  Kravchenko  and  P.  Moin.  Numerical  studies  of  flow  over  a  circular  cylinder  at  Ren  = 
3900.  Phys.  Fluids,  12:403-417,  2000. 

[24]  A.G.  Kravchenko,  P.  Moin,  and  R.  Moser.  Zonal  embedded  grids  for  numerical  simulations  of 
wall-bounded  turbulent  flows.  J.  Comp.  Phys.,  127:412-423,  1996. 

[25]  D.K.  Lilly.  A  proposed  modification  of  the  germano  subgrid  scale  closure  method.  Phys.  Fluids 
A,  4:633-635,  1992. 

[26]  A.L.  Marsden,  M.  Wang,  B.  Mohammadi,  and  P.  Moin.  Shape  optimization  for  trailing- 
edge  noise  control.  Workshop  on  Geometry,  Dynamics  and  Mechanics  in  Honour  of  the  60th 
Birthday  of  J.E.  Marsden.  Toronto,  Canada,  August  7,  2002. 

[27]  P.J.  Mason  and  N.S.  Callen.  On  the  magnitude  of  the  subgrid-scale  eddy  coefficient  in  large- 
eddy  simulations  of  turbulent  channel  flow.  J.  Fluid  Mech.,  162:439-462,  1986. 

[28]  R.  Mittal  and  P.  Moin.  Suitability  of  upwind-biased  finite  difference  schemes  for  large-eddy 
simulation  of  turbulence  flows.  AIAA  J.,  35:1415-1417,  1997. 

[29]  B.  Mohammadi.  Dynamical  approaches  and  incomplete  gradients  for  shape  optimization. 
AIAA  Paper  99-3371  1999. 


46 


[30]  B.  Mohammadi,  J.I.  Moldho,  and  J.  Santiago.  Design  of  minimal  dispersion  fluidic  channels 
in  a  cad-free  framework.  Proc.  Summer  Program ,  Center  for  Turbulence  Research,  NASA 
Ames/Stanford  Univ.,  pp.  49-62,  2000. 

[31]  B.  Mohammadi  and  O.  Pironneau.  Applied  Shape  Optimization  for  Fluids.  Oxford  University 
Press,  Oxford,  2001. 

[32]  F.  Nicoud  and  J.S.  Baggett.  On  the  use  of  optimal  control  theory  for  deriving  wall  models  for 
LES.  Annual  Research  Briefs ,  Center  for  Turbulence  Research,  NASA  Ames/Stanford  Univ., 
pp.  329-341,  1999. 

[33]  F.  Nicoud,  J.S.  Baggett,  P.  Moin,  and  W.  Cabot.  Large  eddy  simulation  wall-modeling  based 
on  suboptimal  control  theory  and  linear  stochastic  estimation.  Phys .  Fluids ,  13:2968-2984, 
2001. 

[34]  A.  Papoulis.  Probability,  Random  Variables,  and  Stochastic  Processes.  McGraw-Hill,  New 
York,  1965. 

[35]  U.  Piomelli  and  E.  Balaras.  Wall-layer  models  for  large-eddy  simulations.  Ann.  Rev.  Fluid 
Mech .,  34:349-374,  2002. 

[36]  U.  Piomelli,  J.  Ferziger,  P.  Moin,  and  J.  Kim.  New  approximate  boundary  conditions  for  large 
eddy  simulations  of  wall  bounded  flows.  Phys.  Fluids  A ,  1:1061-1068,  1989. 

[37]  U.  Schumann.  Subgrid-scale  model  for  finite  difference  simulations  of  turbulent  flows  in  plane 
channels  and  annuli.  J.  Comp.  Phys.,  18:376-404,  1975. 

[38]  W.C.L.  Shih,  C.  Wang,  D.  Coles,  and  A.  Roshko.  Experiments  on  flow  past  rough  circular 
cylinders  at  large  Reynolds  numbers.  J.  Wind  Eng.  and  Industrial  Aerodynamics ,  49:351-368, 
1993. 

[39]  A.  Travin,  M.  Shur,  M.  Strelets,  and  P.  Spalart.  Detached-eddy  simulations  past  a  circular 
cylinder.  Flow  Turb.  Combust ,  63:269-291,  1999. 

[40]  M.  Vainberg.  Variational  Methods  for  the  Study  of  Nonlinear  Operators.  Holden  Day,  San 
Francisco,  1964. 

[41]  M.  Wang  and  P.  Moin.  Computation  of  trailing-edge  flow  and  noise  using  large-eddy  simulation. 
AIAA  J.,  38:2201-2209,  2000. 

[42]  M.  Wang  and  P.  Moin.  Dynamic  wall  modeling  for  large-eddy  simulation  of  complex  turbulent 
flows.  Phys.  Fluids ,  14:2043-2051,  2002. 

[43]  K.A.  Warschauer  and  J.A.  Leene.  Experiments  on  mean  and  fluctuating  pressures  of  circular 
cylinders  at  cross  flow  at  very  high  Reynolds  numbers.  In  Proc.  Int.  Conf.  on  Wind  Effects  on 
Buildings  and  Structures ,  Saikon,  Tokyo,  pp.  305-315,  1971  (see  also  [45]). 

[44]  H.  Werner  and  H.  Wengle.  Large  eddy  simulation  of  turbulent  flow  over  and  around  a  cube  in 
a  plane  channel.  In  Selected  Papers  from  the  8th  Symposium  on  Turbulent  Shear  Flows ,  edited 
by  F.  Durst,  R.  Friedrich,  B.E.  Launder,  U.  Schumann,  and  J.H.  Whitelaw,  Springer,  New 
York,  pp.  155-168,  1993. 

[45]  M.M.  Zdravkovich.  Flow  Around  Circular  Cylinders.  Vol.  1:  Fundamentals .  Oxford  University 
Press,  Chap.  6.,  Oxford,  1997. 


47 


